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Abstract 


We introduce improved numerical techniques for simulating the structural and com- 
positional evolution of planetary rings due to micrometeoroid bombardment and sub- 
sequent ballistic transport of impact ejecta. Our current, robust code is capable of 
modeling structural changes and pollution transport simultaneously over long times 
on both local and global scales. In this paper, we describe the methodology based on 
the original structural code of Durisen et al. (1989, Icarus 80, 136-166) and on the 
pollution transport code of Cuzzi and Estrada (1998, Icarus 132, 1-35). We provide 
demonstrative simulations to compare with, and extend upon previous work, as well 
as examples of how ballistic transport can maintain the observed structure in Saturn’s 
rings using available Cassini occultation optical depth data. In particular, we explicitly 
verify the claim that the inner B (and presumably A) ring edge can be maintained over 
long periods of time due to an ejecta distribution that is heavily biased in the prograde 
direction through a balance between the sharpening effects of ballistic transport and 
the broadening effects of viscosity. We also see that a “ramp” -like feature forms over 
time just inside that edge. However, it does not remain linear for the duration of the 
runs presented here unless a less steep ejecta velocity distribution is adopted. We also 
model the C ring plateaus and find that their outer edges can be maintained at their 
observed sharpness for long periods due to ballistic transport. We hypothesize that 
the addition of a significant component of a retrograde-biased ejecta distribution may 
help explain the linearity of the ramp and is probably essential for maintaining the 
sharpness of C ring plateau inner edges. This component would arise for the subset 
of micrometeoroid impacts which are destructive rather than merely cratering. Such a 
distribution will be introduced in future work. 


1 Introduction 

The rings’ huge surface area-to-mass ratio ensures that they are particularly susceptible to 
the effects of extrinsic meteoroid bombardment. Until recently, the mass of Saturn’s rings 
was thought to be comparable to a Mimas mass (although see Charnoz et al., 2009). Based 
on this value, the area-to-mass ratio for Saturn’s rings exceeds that of Mimas 1 by > 10 4 . A 
consequence of this is that hypervelocity micrometeoroid impacts on the rings, depending 
on parameters, would likely erode them on timescales much shorter than the age of the 
solar system, if all the ejecta escaped (Durisen et al, 1989). More realistically, these impacts 
produce a large amount of particulate ejecta, the vast majority of which are ejected at speeds 
much less than the velocity needed to escape the rings. As a result, a copious exchange of 
ejecta between different ring regions occurs which leads to changes in ring structure and 
composition on both local and global scales (Durisen et al, 1989; Cuzzi and Estrada, 1998; 
Charnoz et al, 2009). This process, by which the rings evolve subsequent to meteoroid 
bombardment, is referred to as “ballistic transport” of impact ejecta (Ip, 1983; Durisen, 
1984a, b; Lissauer, 1984). 

In a series of papers, Durisen and colleagues (Durisen et al., 1989; 1992; 1996) developed 

1 For purpose of this illustration, all the mass is assumed to be in the B ring where the optical depth 
r > 1. The mass of Mimas is ~ 4 x 10 22 g. 
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the first rigorous dynamical code to model ring structural evolution clue to meteoroid bom- 
bardment and ballistic transport. They found that the influence of these processes on the 
rings could explain certain aspects of ring structure such as the fairly abrupt inner edges of 
the A and B rings, and the very similar “ramp” features which connect them to the Cassini 
division and C ring (see Fig. 1) respectively. These structures imply evolutionary times 
of > 100 “gross erosion” times, where the gross erosion time Iq (see Table I for a list of 
parameters) is defined as the time a reference ring annulus of surface density a would dis- 
appear due to loss of ejected material if nothing returned (see below, and Sec. 2.2.1). In 
a complementary study, Cuzzi and Estrada (1998, hereafter CE98) developed a model for 
the evolution of composition while assuming constant structure. They calculated how the 
abundance of both intrinsic and extrinsic non-icy materials evolves over time as icy rings 
are bombarded by largely cometary material, and how these impurities are redistributed 
over the rings. CE98 found that they could simultaneously explain the albedo and color 
dichotomy between the C ring/Cassini division material versus the A ring/B ring material 
and the radial variation of color across the C ring/B ring transition in a time scale similar 
to that on which Durisen and colleagues explained structural evolution. 

Two key quantities in ballistic transport are the impact yield Y and the impacting mi- 
crometeoroid flux <Tj m . Both Y and cr im (see Table 1 and Sec. 2) are essential for providing 
more accurate age-dating of specific ring features, as well as the overall age of the rings 
themselves. The yield of a single impact, Y, which is defined as the ratio of ejecta mass to 
impactor mass, can be quite large depending on several factors (e.g., see Durisen, 1984b). 
The gross erosion time is expressed in terms of these quantities tc = cr/Ya im ; this definition 
serves as a handy reference, but in fact most ejecta are not lost (see Sec. 2.2.1). Moreover, 
this definition is appropriate for cratering impacts, but does not straightforwardly generalize 
to disruptive impacts. 

The ejecta velocity distribution resulting from an impact depends on the hardness of the 
target and the angle of impact (Cuzzi and Durisen, 1990, hereafter CD90). If the target 
is powdery, yields can be in excess of ~ 10 5 — 10 6 for cratering (non-disruptive) impacts 
at normal incidence (e.g., Burns et ai, 1984), while micrometeor-sized particles impacting 
into hard or granular surfaces can have yields as small as rs_/ 10 3 (Vedder, 1972). The ejecta 
velocities for the bulk material from cratering impacts tend to range from ~ 1 — 10 m s _1 , 
much less than the local orbital velocity within the rings (tens of kilometers per second). 
This means that, regardless of whether impacts are cratering or disruptive, the net mass 
gain or loss from the rings due to micrometeoroid bombardment is small compared to its 
redistribution from place to place. The net gain or loss needs to be considered only for very 
long exposure times (Charnoz et al, 2009) or regions where tiny charged ejecta can be swept 
into the planet (Northrop and Connerney, 1987). 

Earlier estimates of the current micrometeoroid flux at Saturn vary significantly (e.g., 
Morhll et al, 1983; CE98; Landgraf et al, 2000), but all suggest that the rings could be 
impacted by close to their own mass (for the Mimas mass estimate) over the age of the 
Solar System. These estimates are largely based on the meteoroid mass fluxes measured by 
the Pioneer and Ulysses spacecraft between 5 — 10 AU (see, CE98, Fig. 17). Some hope 
for improving this estimate had recently surfaced from Galileo measurements of the flux 
at Jupiter. Sremcevfc et al. (2005) used an indirect technique to provide an estimate of 
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the mass flux that may be off by at most a factor of 2 — 3 (larger) compared to previous 
estimates and is comparable to that estimated by CE98 (see Sec. 2.2.1). However, some 
recent measurements of the mass flux at Saturn suggest that the flux may be considerably 
lower than previously thought (Kempf, private comm.), making to longer than what we 
generally assume in this work. 

In a frame rotating at orbital velocity, the ejecta mass from an extrinsic micrometeoroid 
impact that is not disruptive is thrown predominantly in the prograde orbital direction. This 
result arises naturally from the tendency of projectiles to arrive with velocity vectors that 
are retrograde in the orbital sense, as azimuthally averaged over the rings (CD90). That is, 
the combination of the orbital motion of ring particles and the motion of Saturn through 
the micrometeoroid flux leads to a “headwind” of material that increases both the impacts 
on the leading faces of ring particles and the impacting velocities (see also Latter et al, 
2012). In addition to mass, ejecta carry away with them angular momentum. Since most 
of the ejecta from a non-disruptive (cratering) impact are prograde, they tend to reimpact 
the rings at outer locations. Prograde ejecta are launched from their original radial location 
with more angular momentum than their parent ring particle, but less angular momentum 
than ring particles they may impact on their next ring crossing at some outer radius. The 
net result is to decrease angular momentum at the secondary impact location, leading to 
radial inward drift. 

On the other hand, an impact that leads to complete disruption of a target ring par- 
ticle into several fragments would likely produce the opposite effect because the fragment 
velocities are biased in the same direction as that of the impactor, resulting in a retrograde 
distribution (relative to Keplcrian) with lower ejecta velocities than their prograde counter- 
parts (Nakamura and Fujiwara, 1991; Paolicchi et a/., 1989). In either case, the structure of 
the rings (he., optical depth r and surface density a) will have an effect on the rate of mate- 
rial drift. This is because the probability of ejecta absorption, which determines the actual 
ejecta mass that is reabsorbed by the rings instead of merely passing through them, depends 
on the local r, and its angular momentum depends linearly on a (CD90; CE98; Charnoz et 
al, 2009). In fact, there is a smaller, but systematic, tendency for infall to decrease the net 
angular momentum of rings (Ip, 1983; CD90; see below) 

These processes have been used to estimate ring lifetimes in several ways. First, they 
will lead to angular momentum loss within the rings globally. For example, it has been 
estimated that the C ring would be lost to the planet in rs_/ 10' — 10 8 years (CD90) if no other 
mechanism were at work to sustain it. A similar age for the rings related to the loss of grains 
and molecules to Saturn’s ionosphere generated by impacts was obtained by Northrop and 
Connerney (1987) 2 . Second, carbon and silicate-rich meteoroid material also darkens and 
“pollutes” the almost entirely icy rings over time, an effect that provides a powerful means 
for estimating ring age. Doyle et al. (1989), and subsequently CE98, using a mass flux from 
Grim et al. (1985), noted that the relatively high albedo of the A and B rings was inconsistent 
with their having retained more than a small fraction of primitive, carbonaceous material 
they would have accreted over the age of the Solar System as compared to the expected 
fraction of order unity. CE98 further demonstrated using their pollution transport model 

2 It should be noted that resonant interactions with nearby ring- moons can also lead to significant angular 
momentum loss (e.g., Goldreich and Tremaine, 1982; Esposito, 1987; Poulet and Sicardy, 2001). 
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that the relatively small amount of extrinsic darkening material needed to evolve the inner 
B ring and the C ring to their current broadband colors would also suggest a geologically 
young age for the C ring, similar to the time it would take to lose the ring based on angular 
momentum loss arguments. This would seem to pose a problem for the possibility that the 
rings are more massive than previously thought, unless the C ring we see today is some 
great-grandchild of an original ancestor (Charnoz et ai, 2009, see Sec. 3.1). 

In fact the C ring provides fertile ground for ballistic transport modeling because most 
of the structure in the C ring remains unexplained, and there are several candidates for 
ballistic transport signatures. For instance (see upper panel of Figure 1), Durisen et al. 
(1992) suggested that the sharp inner edge of the B ring was maintained by a balance 
between ballistic transport, which tends to sharpen low-to-high optical depth transitions, 
and the broadening effect of viscous transport. Likewise, these workers also found, using a 
prograde ejecta distribution (CD90), that the adjoining “ramp” structure on the low optical 
depth side of the B ring inner edge inevitably came about due to advective effects. A similar 
sharp edge connected to a ramp is present at the inner A ring/Cassini division boundary 
(Fig. 1, lower panel). In this paper, we examine these results in more detail and establish 
that these features are indeed products of ballistic transport. 

The shapes of the C ring “plateau” peaks are also suggestive (see Colwell et ai, 2009). 
Material tends to drift more quickly in low optical depth, featureless regions because the 
kinematic viscosity is weak there. This means that inwardly drifting material will begin to 
“pile up” which leads to a sharpening of an outer plateau edge. Likewise, retrograde ejected 
material would exhibit outward radial drift and would tend to lead to a pile up at the 
plateau’s inner edge. A plateau might then exhibit a characteristic in which an outer (inner) 
edge has higher amplitude than the inner (outer) edge, with a ramp in optical depth between 
them, if prograde (retrograde) ejecta is the dominant component. The observed plateaus in 
Fig. 1 display these features to differing degrees. Interestingly, there appears to be some 
symmetry between inner/outer plateau edge structure about the Maxwell ringlet at around 
87500 km (~ 1.45 Rs, where Rs = 60330 km is the radius of Saturn). Whether this is a 
coincidence or not is unclear. It may be that some variable combination of particle size, which 
may even differ between plateaus and surrounding regions (Colwell et ai, 2012; Marouf et 
al, 2012), and impact velocities may lead to a transition between primarily cratering impacts 
(prograde ejecta) versus primarily disruptive impacts (retrograde ejecta) at around 1.45 Rs 
(Estrada and Durisen, 2011). 

Despite these intriguing possibilities, the plateaus, as well as the undulating structure 
in the inner C ring and the semi-regular ~ 80 km structure in the inner B ring, none of 
which are associated with known resonances (although see Hedman and Nicholson, [2013]), 
have yet to be modeled successfully in convincing detail. Although we do not study the 
issue explicitly in this paper, an instability associated with ballistic transport (hereafter, 
BTI, Durisen 1995; Latter et a/., 2012; 2014a, b) can lead to the spontaneous creation of 
undulatory structure in an otherwise uniform ring if ring conditions fall within a specific 
range of optical depths and wavelengths. This process could play an important role in the C 
ring and inner B ring. Specifically, Durisen (1995) suggested that ballistic transport could 
cause the semi-regular structure seen in the inner B ring and possibly the longer wavelength 
C ring undulations without the inner edge being the driving mechanism. Latter et al. (2012; 
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2014a, b) were able to confirm and extend on the analysis of Dnrisen (1995). The simulations 
we present in this paper do strongly suggest that the BTI is at work, but operating under a 
very nonuniform background structure, where it is not as well studied or understood. 

Another important issue, which we briefly touch upon in this paper, is whether the rings 
of Saturn are much more massive than previously thought (Sec. 3.2). If so, high optical 
depth regions, i.e., the B ring where most of this extra mass must be hidden (e.g., Charnoz 
et al, 2009), should be able to resist pollution from extrinsic bombarding material even more 
easily than previously thought, and this should be manifested in current particle albedos. 
A more massive B ring does not necessarily indicate that all the ring elements are old, and 
indeed, the C ring probably cannot be long-lived. Charnoz et al. (2009) argue that the most 
likely time that the main rings formed would have been at the Late Heavy Bombardment 
(LHB, Gomes et al, 2005), setting up a mismatch between the age of the main rings and 
the age of the C ring. 

A review of previous work on the subject of ballistic transport can be found in Charnoz 
et al. (2009). In this paper, we lay the groundwork for more detailed future modeling of 
micrometeoroid bombardment and ballistic transport in planetary rings through the intro- 
duction of an improved ballistic transport code that is capable of modeling both structure 
and composition in tandem. With this code, we are able to model the rings over a broad 
range of spatial scales and over both short and long timescales. These are important devel- 
opments because the results of this work in combination with that of Latter et al. (2012; 
20 14a, b) strongly suggest that ballistic transport is operating within the rings, but it will be 
a more complicated endeavor to match all the observed features than previously envisioned, 
requiring additional physics that we have not yet incorporated in our model. Our improved 
code, though, is robust enough that new physics can readily be included, which will afford 
us the ability to explore this complexity in more detail. 

In Section 2, we describe the development of our combined structural and compositional 
code. In Section 3, we present the results of demonstrative simulations which reproduce and 
extend upon previous work, as well as introduce new modeling of observed C ring structure. 
In Section 4 we present a discussion of our conclusions and, finally, in Section 5 we outline 
the direction of future work. 


2 Combined Structural and Compositional Model 

2.1 Limitations of Previous Work 

Previous studies involving ballistic transport in Saturn’s main rings were limited in two ways. 
From the structural standpoint, solving the ballistic transport equations (see Sec. 2.2.3) is 
computationally expensive, even with today’s processing abilities. Ballistic transport essen- 
tially requires tracking the exchange of material from one location (radial bin) in the radial 
range of computation to every other radial location within its maximum “throw distance” , 
which depends on the upper limit of the ejecta velocity distribution. Meanwhile, one must 
also sum up all the contributions to that radial bin from all other radial locations from which 
ejecta can reach it. Because this involves integration over multiple variables (e.g., angles, 
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velocities, reabsorption probabilities) at every time over a potentially large histogram of ra- 
dial bins, comprehensive parameter studies were not feasible, and global scale calculations 
intractable. 

Another important handicap in previous studies was a lack of adequate data for some 
basic ring properties. For example, prior to Cassini, the ring optical depth t(R) was known 
in most locations, but the ring surface density <j(R) had only been measured in several 
isolated spots, which were generally associated with spiral density wave regions. It was then 
only possible to determine the ring opacity k(R) = r(R)/a(R) in these discrete locations 
(e.g., Cuzzi et al, 1984; Lissauer and Cuzzi, 1985). CE98 extended these values to non- 
density wave regions, assuming wave regions were not too different from non-wave regions, 
but could only estimate the spatial variability of k in optically thick regions because of a 
lack of occultation data (see Sec. 2.2.6). All of the simulations of Durisen et al. adopted 
the condition r oc a and thus the opacity was assumed to have a constant value dependent 
on some reference value of a in the inner B ring (Sec. 3.1, also see Cuzzi et al, 1984). 

A major limitation of the C ring/inner B ring compositional study of CE98 was their 
assumption that the underlying ring structure could be kept constant while the ring com- 
position evolved, based on the finding by Durisen et al. (1992) that constant optical depth 
regions and inner sharp edges can remain more or less unchanged over long timescales. As a 
result, CE98 assumed that the optical depth and mass surface density were time-independent 
on the timescale of compositional evolution. This meant that the cumbersome r-dependent 
loss and gain integrals (see Sec. 2.2.5) that determine the redistribution of mass and an- 
gular momentum due to ballistic transport needed to be done only once which provided a 
considerable savings in computational time. 

Cassini observations are providing a more detailed understanding of dynamical ring prop- 
erties such as the opacity k, surface density a, and optical depth r which are key quantities 
in modeling of ring structure. These observations, combined with numerical improvements 
will allow us to embark on much more detailed modeling and analyses than could have ever 
been done previously. On the other hand, Cassini has also shown us how local measurements 
of r suggest a bimodal or multimodal “gap and clump” structure (he., self-gravity wakes) 
that we do not deal with explicitly here (e.g., Colwell et al . , 2009) but hope to incorporate 
into future simulations. 

2.2 The Ballistic Transport Code 

The ballistic transport code (hereafter BT code) evolves the surface mass density of the rings 
a(R,t), and thus its structure, over time by solving the system of equations given below in 
Sec. 2.2.3. The code also tracks the changes in the fractional mass of extrinsic (bombarding) 
material, or “pollutant” , which is assumed to be made up of some fraction of icy and non-icy 
constituents, and the rings’ intrinsic material which is assumed to be composed mostly of 
water ice mixed with a very small fraction of absorbing material which is different from 
the micrometeoroid impactors. Thus the BT code also simultaneously evolves the bulk 
composition of the rings. 

In ballistic transport, the two main effects that contribute to the rings’ evolution are the 
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direct net mass exchange of micrometeoroid impact ejecta between different ring regions (Ip, 
1983), and the “indirect” compression or expansion arising from differential radial drift of 
material due to the direct net angular momentum exchange associated with the redistribution 
of impact ejecta (Lissauer, 1984; CD90). Secondary effects (e.g., mass loading, ring torque) 
that are related to the direct deposition of extrinsic material and its redistribution can also 
play a role but are lesser in magnitude and more relevant over long timescales (Durisen et 
al, 1996; see Sec. 2.2.4). 

The exchanges in mass and angular momentum h per unit mass are determined at each 
radial location R in a Lagrangian radial grid (he., one that moves with the material) by 
computing (1) the local gain rate in mass and angular momentum attributed to all of the 
extrinsic material and incoming ejecta from neighboring ring regions that are absorbed at R 
and (2) the local loss rate in both quantities by accounting for the all of the ejecta created 
at R , which are ejected at different speeds and in different directions, that are absorbed 
at radial locations different from R. The net changes in o and h at a given R are then 
determined by differencing the local gains and losses. We call this combination of effects the 
“direct” contribution. Advection further modifies the local surface mass density through the 
compression or expansion of the Lagrangian cells based on differential radial drift velocities 
due to angular momentum readjustment, viscous angular momentum transport, and several 
other sources of radial drift. We call this the “indirect” contribution. 

In this section we describe the motivation and mathematical formulation for the BT 
code. The most relevant parameters used are listed in Table 1. Much of what is presented 
below has been described elsewhere in detail (e.g., Durisen et al, 1989; 1992; 1996; CE98); 
however, given the mathematical and numerical focus of this paper, we believe it is justified 
to repeat some of these previous descriptions, while emphasizing differences in approach and 
new physics when appropriate. More detailed discussions of some of the various aspects of 
the BT code are included in the appendices; in particular, the numerical methodology of the 
code itself is described in Appendix D. 

2.2.1 The impact flux on the rings and definition of the gross erosion time 

Of the various input parameters for the modeling of ballistic transport in Saturn’s rings, the 
micrometeoroid flux remains the most fundamental. For this work, we continue to use the 
value of the one-sided, incident flux at infinity = 4.5 x 10 17 g cm 2 s 1 used in CE98. 
This value was initially obtained from direct integration over the 1 AU interplanetary flux 
models of Griin et al. (1985) by CD90, and was smaller than the value used by either Morfill 
et al. (1983) or Ip (1984). This value was used subsequently by Durisen and colleagues in 
later papers. CE98 reanalyzed and corrected the estimates from CD90 by including Pioneer 
10/11 and Ulysses spacecraft data (see Fig. 17, CE98). CE98 found that the CD90 value 
for boo remained plausible within a factor of 2 — 3, in spite of an inconsistency in the CD90 
assumptions (see CE98 for a detailed discussion). The two-sided flux of micrometeorites 
impacting the ring plane and potentially absorbed at some radial location R is given by &i m : 

& im = 2 & oo ^(t)F g (R/R 0 )-°- s , (1) 
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where we have accounted for the gravitational focusing of the planet (see CD90) using an 
asymptotic relative velocity of 14 km s -1 , which is numerically averaged into a factor F g = 3 
at a reference radius Ro = 1.8 R5. The flux at other radii is obtained by the numerical fit 
(R/Ro)~ 0 ' 8 to the radial dependence of the calculated focusing (CD90; also see CE98). The 
optical depth dependent function 

s# (r) = [1 - e- (T/Ts)P ] Vp , ( 2) 

is a parameterization of a numerically determined impact probability that depends on r and 
the various angular aberrations averaged over the orbit of the ring particles (CD90; Durisen 
et al, 1996). Here the fit parameters are t s = 0.515 and p = 1.0335 (CE98). 

The actual mass flux of ejected material <j e j at reference radius i? 0 is determined by 
integration over the full ejecta yield function <3/ derived by CD90 (see their equations 40- 
46) 3 : 


cr ej (i?o) = 2(700 J J 3^(R 0 ,x, tydxdtt « Ya im (R^), (3) 

which depends on the ratio of the ejecta velocity to orbital velocity x = v^/vk, and solid angle 
Q(9, (j)), where 9 and (f) are the zenith and azimuth angles of emergent ejecta as measured from 
the ring particle velocity vector in the frame of the planet and the ring normal. The ejecta 
distribution function 3/ contains all the information about impact and escape probabilities 
and about gravitational focusing and is normalized to a yield of Y = 10 4 at 14 km s _1 
impact velocity (Lange and Ahrens, 1987). We use the approximation on the RHS of Eq. 
(3) below for illustrative purposes. It typically differs by no more than tens of percent from 
the exact calculation obtained through direct integration of , (CE98, see Appendix A) 
which we always utilize in the BT code. 

The “gross erosion time” t G is the fundamental time unit of ballistic transport and is 
defined as the time it would take for a ringlet of surface density a to be completely eroded 
away if no material returned. 


tG 


0 a 





( 4 ) 


For example, using the two-sided, gravitationally focused mass flux in Eq. (1) and an ejecta 
yield of Y = 10 4 , a reference ring annulus with a = 100 g cm -2 would erode away in t G > 10 6 
years. 

We normalize the time units of the BT code to the gross erosion time at a particular 
radius, allowing our simulations to be independent of specific values of the micrometeoroid 
flux and yield. For specificity and ease of comparison, the evolutionary times for this paper 
are in the same units as those used in past papers (e.g., Durisen et al. , 1992; CE98) 

3 Note that in CD90, the fitted solution to 'W is expressed in terms of the cone and clock angles of the 

emergent ejecta. 



years, 


( 5 ) 


t G ~ 1.3 x 10 5 



4.5 x 10 17 g cm 2 s 1 


(Ti, 


°( T 

96 g cm -2 ) 


which has been defined for r = 1 at a radius of R — 1.52 R5 which corresponds to the B 
ring inner edge. Defining t G at this radius was sensible for the edge simulations of Durisen 
and colleagues, but one must keep in mind that the gross erosion times for the lower optical 
regions of the C ring (which range from r ~ 0.05 — 0.5) will be shorter than the time given in 
Eq. (5). In any case, the absolute ages derived from our simulations remain uncertain because 
of uncertainties in &i m , Y, and the efficiency with which extrinsic bombarding material retains 
its absorptive properties. 


2.2.2 Code assumptions 

Basics. The rings are modeled as a Keplerian thin disk that is axially symmetric, with the 
central planet taken to be a point mass, and ring self-gravity is ignored. A thin disk is a 
reasonable assumption because the velocity dispersion c of ring particles is typically a few 
millimeters per second (Schmidt et al, 2009; Colwell et al, 2009) somewhat akin to the 
escape velocity from several- meter- sized particles, which are the sizes observed to be the 
upper bound in the particle size distribution (Cuzzi et al, 2009). On the other hand, the 
orbital speeds of ring particles Vk are on the order of tens of kilometers per second. The 
random motions clue to c lead to a vertical thickness H/R ~ c/% C 1, or an H of a few 
times our largest particle size (Cuzzi et al, 1979). 

To consider the disk to be thin with respect to ballistic transport, the radial and vertical 
excursions of impact ejecta should be much greater than H . Ejecta velocities n e j from non- 
disruptive impacts can vary from ~ 1 — 100 m s _1 depending on the hardness of the target 
(e.g., Lange and Ahrens, 1982; Hartmann, 1985). The inward or outward radial excursion of 
an ejecta particle, which we refer to as its ‘throw distance’, is given by \SR\ « 4 xR (Morfill 
et al, 1983; Durisen et al, 1992). Comparing the radial throw distance to the ring vertical 
thickness, \SR\/H ~ 4 xR/ R(c/vk) ~ 4u e j /c 3> 1. By the same token, the average vertical 
excursion is \8z\ ~ v^P/2 ~ ttR^v^/vk), where P is the orbital period, so that again we find 
\8z\/H ~ vruej/c S> 1. 

Our assumption of axial symmetry stems from a comparison of timescales. The typical 
orbital period of a ring particle is on the order of ten hours, whereas the relevant timescales 
of ballistic transport are on the order of t G ~ 10 4 — 10 6 years. Given this disparity, it 
is reasonable to expect that transient deviations from axial symmetry in underlying ring 
structure due to impacts average out over the timescales involved in our simulations. We 
use a variable timestep in our BT code (see Appendix D) which could be as short as a few 
years if the ring ambient conditions in the simulation call for it. This still equates to > 10 3 
orbits which should not invalidate our assumption. We leave consideration of the much 
more complicated problem of 2-D ballistic transport modeling (e.g., allowing explicitly for 
non-axisymmetric wakes [Colwell et al, 2009]) for future study. 

Ring self-gravity is considered to be unimportant with respect to impact ejecta because 
the escape velocity from the ring particle size that dominates the Saturn ring mass, a few 
meters, is much smaller than the typical ejection speeds we use here. The presence of larger 
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moonlets whose masses are large enough to produce escape velocities on the same order as 
that of the typical ejectum are also not considered, because these objects are sparse and thus 
would only effect a very small fraction of ejecta at any given time. Saturn’s non-sphericity 
can lead to nodal precession and apsidal motion on timescales of tens of days in the C 
ring, with timescales increasing with radial location (Lissauer and Cuzzi, 1982). Given the 
relatively short orbital times of ring particles, ejecta are expected to complete at most a 
few orbits before being reabsorbed (Durisen et al, 1989; CD90). Thus, Saturn’s oblateness 
should make little difference for ejecta trajectories and is neglected. 

Ejecta distribution. We assume that all ejecta travel on Keplerian elliptical orbits between 
their point of ejection and their point of eventual reabsorption, and reabsorbed ejecta do not 
produce secondary ejecta distributions. The former assumption is valid because x <C 1, so 
the ejecta travel on slightly elliptical orbits and cannot escape the system but rather reimpact 
the ring at distances up to the maximum throw distance ~ 4 x max R and in timescales much 
less than the orbit precession timescale. The latter assertion seems sensible in that the 
secondary impacts of ejecta grains with the grainy regoliths of ring particles would be at less 
than ~ 100 m s _1 velocities instead of tens of kilometers per second and likely to be highly 
inelastic (e.g., Hartmann, 1985). The secondary ejecta velocities are thus expected to be 
•C v e j, with reimpacts occurring locally, within the maximum radial bin resolution we use in 
the BT code (> a few kilometers). 

Ring and meteoroid composition. The intrinsic composition of the rings is assumed to be 
primarily water-ice with a small fraction of spectrally red, absorbing material that is “intra- 
mixed” within the ice grains. CE98 found that Titan tholins provided a good fit to the 
observations, however other candidates for the reddish “UV-absorber” have been proposed 
(e.g., Cuzzi et al., 2009; Clark et al, 2012). The extrinsic, bombarding material is assumed 
to contain a fraction / ext ~ 0.5 of a spectrally neutral, non-icy darkening agent, with the 
remainder being water ice. The BT code tracks the mass fraction f e (R, t ) of non-icy extrinsic 
material in the rings and assumes that the remainder (1 — / e ) is the admixture of primordial 
and extrinsic icy component. Extrinsic absorbing material is assumed to be retained locally 
in the ring with an efficiency r) (see, e.g., Doyle et al., 1989; CE98), and uniformly mixed 
throughout the local mass density cr(R,t). Thus, all the components of composition are 
assumed to be volumetrically mixed within the ring particles at all times. 

This assumption is more than likely adequate over long timescales (> a few tc)- However, 
a more recent model of regolith growth in ring particles due to bombardment by a broader 
size range of projectiles (up to meter size) suggests that, for ring particles that initially begin 
as pure ice, their regoliths remain more or less constant in the fractional amount of pollutant 
for timescales < 1 tc as they grow in depth (Elliott and Esposito, 2011). Once the regolith 
has fully developed and essentially fills the particle, fractional mass increases linearly with 
time which is consistent with this work. The work of Elliott and Esposito (2011) becomes 
particularly relevant when considering the study of shorter timescale, transient features. For 
example, a recent local event such as the breakup of a small icy moon may lead to the 
resetting, locally, of the ballistic transport “clock”, an idea that has already been suggested 
to explain variations in brightness in the rings over length scales of ~ 1000 — 3000 km that do 
not seem consistent with ballistic transport acting over long timescales (Cuzzi et al., 1984; 
Esposito et al, 2005). In future work, we plan to incorporate these effects into our BT code. 
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In general though, the gardening process modeled by Elliott and Esposito (2011) validates 
our assumption that a ring particle’s composition is well mixed and not characterized by 
their surface veneers. 

Conservation of energy. Finally, the ballistic transport equations (Sec. 2.2.3) account 
for conservation of mass and angular momentum, but not of energy (Durisen et al, 1989). 
The reason is that micrometeoroid impacts play little if any role in maintaining the velocity 
dispersion in the rings against damping due to collisions, on energetics grounds. To see this, 
we consider the impact of a projectile of mass m ; m with a ring particle of mass m. Through 
momentum exchange arguments, the random velocity imparted to a ring particle m through 
the impact with a particle rn im is Av ~ Y /i m w e j , where /i m is the mass ratio m im /(m + m im ). 
The ring particle m is impacted on a timescale 


fim « « thn — = (/im Y) t G . (6) 

di m 

Then, the rate at which specihc energy is introduced into the system as a result of these 
impacts is Av 2 /ti m , which can be compared with the rate at which specihc energy is dissipated 
due to inelastic collisions (Cuzzi et al, 1979) 


Av 2 v% c 2 rfl(l - e 2 ) 

— — = ffirA T 1 > ^ 

6im t'G 


hr 


> 1.25 x 10 


_ 12 c 2 ar(l - e 2 ) 


x * 


( 7 ) 


where e is the coefficient of restitution 4 . Using our canonical choice for upper size cutoff 
r 2 = 3 m as the size for m and c ~ 0.2 cm s'” 1 and considering typical parameters in our 
simulations for a C ring plateau, r = 0.4, a — 40 g cm -2 and x = 1CU 3 , one Ends that 
r im > 3(1 — e 2 ) 1 / 3 cm. The mass peak in the interplanetary micrometeoroid flux, which 
is presumed to be cometary in origin, is around a size of 100 /im (e.g., Grim et al., 1985; 
CE98), while Edgeworth-Kuiper belt dust particles have a mass peak between sizes of 1 — 10 
/im (Poppe et al, 2010; Poppe and Horanyi, 2012). Although ring particle collisions are 
moderately elastic, coefficients of restitution very close to unity would be needed for impacts 
by micrometeorites to be able to maintain the observed velocity dispersion in the rings. Thus 
we assume that the velocity dispersion in the rings is maintained by a balance of collisional 
energy loss and viscous gain, which we do not model, and, as a consequence, we do not 
consider conservation of energy in our system of equations (Sec. 2.2.3). 


2.2.3 Working equations 

Given our set of assumptions from Sec. 2.2.2, the ring system is described at any time t 
by a surface mass density a(R,t) and normal optical depth r(R,t), where the two quan- 
tities are related through the ring opacity n(R,t) = r/cr. The opacity may explicitly or 
implicitly incorporate information or assumptions about the ring particle density p and size 
distribution. 

4 Alternatively, one can compare the specific energy introduced into the system to the viscous dissipation 
energy rate (e.g., Araki and Tremaine, 1986; Latter and Ogilvie, 2008). 
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The evolution of the ring system under ballistic transport is determined through solving 
the following set of equations. The mass continuity equation 


da 

~dt 


IA 

RdR 


(aRv R ) = T 


A m T *^ii 


or, alternatively, the areal angular momentum density continuity equation 


( 8 ) 


d_ 

dt 


(ah c ) + 


]__d_ 

RdR 


R(ah c v r) — auR 3 


_d_ 

dR 




(9) 


where h c = Rvk = ( GMR ) 1//2 is the specific angular momentum for circular motion, G is 
the gravitational constant and M is the mass of the planet. Equation (8) follows the local 
changes in a and R due to the net loss or gain of ejecta to and from other ring regions. The 
RHS of Eq. (8) accounts for the direct mass gains T m and losses A m by integrating the ejecta 
distribution over the rings taking into account cylindrical effects, while the RHS of Eq. (9) 
are the angular momentum loss and gain integrals (see Sec. 2.2.5). The divergence of the 
radial mass flux in the ring plane caused by the slow radial drift of ring material at velocity 
Ur brought on by angular momentum transport (indirect contribution) is accounted for in 
the second term on the LHS of (8). 

Changes in the composition of the rings in terms of the mass fraction of non-icy, extrinsic 
material f e (R , t) are followed using the equations for the evolution of the surface mass density 
of the extrinsic and intrinsic components 


da e 

dt 


1 d 

+ — 7^ (a e Rv r) = r/./extbim + r, 


fe(r,t) Am 


Vfext&im T A(7 e , 


( 10 ) 


d(j‘ 1 d 

~dT + ~R~dR ( a i Rv R ) = i 1 - V f ext) dim + - [1 - f e (r, f)]A m = (1 - + A & h (11) 

where a = a e + op Unless otherwise noted, for the simulations in this paper we choose 
/ext = 0.5 and 77 = 0.12 which is consistent with CE98. The radial drift velocity Ur is 
discussed in the next section. The gain integrals T m)e and T mii are described in Sec. 2.2.5. 

2.2.4 The Advective Term 

The total radial drift velocity Ur = Up a11 + u^ lsc + u ^ iasl + u^ rq that appears in the indirect term 
(e.g., second term on the LHS of Eq. [8]) has multiple contributions, the most important 
being that due to the ballistic transport mechanism itself: 

[T h - A h - h c (T m - A m )] . (12) 

Equation (12) can readily be derived from combining Eqns. (8) and (9) for the inviscid 
case where mass loading and net torques are ignored. A second component of radial drift 


12 



is clue to the viscous angular momentum transport. For a thin viscous Keplerian disk 
(Hartmann, 2009), 

= (13) 

where v(R,t) is the kinematic viscosity, which is discussed in detail in Sec. 2.2.7. 

Secondary contributions to the radial drift of material that are due to the collective effects 
of meteoroids (mass loading and asymmetric absorption of ejecta, Durisen et a/., 1996) are 
also included in the calculation of the total Ur in our BT code. Meteoroid bombardment 
leads to the direct deposition of mass, which is accounted for through d im . However a more 
important effect of this bombardment is the change it imparts to the local specific angular 
momentum h. An annulus subject to these effects will drift radially, adjusting its h to match 
that of the equilibrium h = h c at a different radius, with drift velocities due to these collective 
effects Ur 11 given by (Durisen et al., 1996; cf. Eq. [12]) 

„,coll torq , masl ( \ /, V O 

% =»R +% ={<’ 4 r ) [UJ lm _ c1 T (14) 

where the Erst term in the RHS of (14) is the meteoroid angular momentum deposition 
rate per unit area and time, and the second term on the RHS is the mass deposition rate. 
Secondary contributions are discussed in more detail in Appendix B. 

2.2.5 Computation of losses and gains 

The essence of the ballistic transport process lies in the changes of the local surface mass 
density due to the differences between losses and gains of ejecta material at some R and 
the resultant drifts. The differences are of order ~ x 2 away from sharp edges, and ~ x at 
sharp edges, and there are order ~ x 2 differences between the radial velocities wr across 
neighboring ring regions through the divergence term. Quite literally, ballistic transport 
depends on differences of small differences between losses and gains (Durisen et al, 1989). 
The losses A and gains T in mass and angular momentum per unit area and time must be 
computed every time step, and at each radial location within the ring. The losses in mass 
and angular momentum must account for the total rate at which material and its associated 
angular momentum is being carried from ring radius R , in all possible directions and at all 
possible speeds, to other locations R r where the ejecta are reabsorbed: 

POO /*27T p1 

A m = M{R,r) / dx dtp j ^(R, R')F(x,9,(p) d cos9, (15) 

Jo Jo J - 1 

poo p 2n p 1 

Ah = h c M(R, t) / dx d(f) I (1 + xcos9)£P(R, R')F(x, 9, (f>) dcosO. (16) 
Jo Jo J-i 

Here, M gives the total rate of mass ejection at radius R per unit area, and 9 and (j) have been 
defined previously in Sec. 2.2.2. The function F(x, 9, (J) specifies the distribution of ejecta 
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that escape R, as a function of velocity and solid angle defined by 9 and p. By definition, 
the integral fff F(x,9,p) dx dp d cos 9 = 1. The probability function &(R,R') is given by 


&>{R,R') 


1 — e 

1 — e -(r z +T') ’ 


( 17 ) 


where r z and r' are the slant optical depths for the angles of incidence in the rotating frame 
measured from the ring plane normal at R and R', respectively (see Eqns. [29] and [30], 
Durisen et al, 1989). Equation (17) gives the total probability that the ejected particle will 
be absorbed at R' rather than at R. 

The corresponding gains in mass and angular momentum per unit area and time at the 
radial location R are calculated with the integrals T m and Th, which account for all the 
ejecta from other ring regions R! that reach R and are absorbed there. Because ballistic 
transport is sensitive to the small net differences between the losses and gains, Durisen et 
al. (1989) argued that the gain integrals should be as closely analogous to Eqns. (15) and 
(16) as possible for purposes of numerical accuracy. This was accomplished by relating the 
radius of ejection R to all possible radii of absorption R' using a mathematical description 
of the Keplerian elliptical orbit followed by an ejectum R' = RA(x, 9 , p ), where the function 
A is given by Eq. (33) of Durisen et al. (1996). Summing up all of the contributions to the 
ring annulus 2nR dR from all annuli 2ttR 1 dR 1 then leads to the following expressions for the 
mass and angular momentum gain integrals 


r 


m 




3g(K,T')A 2 &’(R!,K)F'(x,e,<l>)d cos 0 


(18) 


POO p2lt pi 

Th — I dx dip h' c {l + x cos 9)&(R', t , )A 2 ^ > (R', R)F\x, 9, p) dcos9, (19) 
Jo Jo J - 1 

where quantities within the integrals are evaluated at R', which is itself evaluated in terms 
of x, 9 and p. The gain integral T m can be broken into the fractional mass gain integrals 
r m e and r m>i , which incorporate the mass fractions of extrinsic and intrinsic materials, but 
are nonetheless similar to Eq. (18). For example, 

poo p2n pi 

r m , e = / dx dp f e {R', t)& (R', F) A 2 A* (R',R)F\x, 9, p)d cos 9, (20) 

Jo Jo J - 1 

and the analogous expression for T mji follows from replacing f e (R',t ) with 1 — f e (R',t). In 
Appendix D we describe our technique for solving these integrals numerically. 

2.2.6 Ring opacity models and particle size 

Knowing both the surface mass density and the optical depth allows one to calculate the ring 
opacity k. In the work of CE98, who investigated whether pollution transport alone could 
explain the detailed shape of the ring color prohles across the rings, only two such locations 
were known in the inner B ring and the C ring: the Janus 2:1 density wave and the Maxwell 
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ringlet, respectively. In order to achieve their goal, CE98 required an opacity profile with 
higher radial resolution. They assumed the particle size distribution is best described by a 
power law in radius r p (R , t) oc r -3 in most regions with upper and lower cutoffs rq, rq of 3 cm 
and 3 nr, respectively (Zebker et al, 1985; Cuzzi et a/., 2009), and we continue to use this as 
our canonical size distribution. Given this distribution, k depends primarily on the largest 
particle size 5 . Under this assumption, CE98 (see their Fig. 5) utilized the variance technique 
of Showalter and Nicholson (1990), which determines the largest “particle” size with high 
radial resolution, to determine from the h-Sco Voyager stellar occultation the largest effective 
radius particle across a range of ring radii where no spiral density waves existed. The scatter 
in their derived hi values, coupled with values inferred by others, led CE98 to adopt very 
simple k(R) profiles for the C ring and inner B ring simulations: 


{ K\ R < Ri, 

( R - R i) Ri<R<R2 (21) 

H 2 R A R-2- 

where (in units of cm 2 g^ 1 ) K\ = 0.05, K 2 = 0.009, Ri = 1.5 Rs (roughly the base of the C 
ring ramp, see Fig. 1), and i ?2 was modestly varied to different locations in the very inner 
B ring. The estimated n = H 2 for the B ring opacity was actually an average over a range 
of possible values, due to a lack of occultation data in optically thick regions (see Fig. 11 of 
Showalter and Nicholson 1990). Recall that the treatment of CE98 (and below) differs from 
that of Durisen and colleagues who effectively used a constant opacity in their models (e.g., 
see Durisen et al, 1992; Sec. 2.1 and 3.1). 

One observation that CE98 gleaned from point-by-point inspection of their derived opac- 
ity values in the C ring was that these values indicated a preferential association of low n 
(larger particles) with the relatively higher optical depth plateaus, while higher k were as- 
sociated with regions away from plateaus where r was lower. This was somewhat supported 
at the time by comparison of Voyager visual and radio wavelength optical depths. If true, 
this effectively would say that there is more mass in the plateaus than the simple model in 
Eq. (21) would predict. This suggests an alternative model for the ring opacity that we 
utilize for this work, namely that the opacity depends on r (or alternatively, a). For models 
that explore the C ring and inner B ring, we can define the initial values of the opacity 
k(R, 0) = a + b/r(R, 0) and at all later times through 


k(R, t ) 


a a 4 b 

2 + 2 [ 1 + a 2 a{R,t)\ ’ 


( 22 ) 


where the fit parameters a — K\ + (A/t/ A t)t 2 , b = — (An/ A t)t\T 2 , Ak = k 2 — n\ and 
At = 72 — T\. The values of and r ± t 2 can readily be chosen to be consistent with the 
profile given in Eq. (21). Thus, for T\ = 0.05 and r 2 = 1, a ~ 0.0068 and b ~ 0.0022 (both 
cm 2 g _1 ), respectively. The above formula is heuristic, but it allows us to associate a certain 
optical depth with a specific surface mass density. This in turn allows us to model plateaus 


5 For power law slopes > 3, this assumption is not true because the lower bound of the particle size 
distribution becomes increasingly important in determining n (e.g., see CE98). 
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that are relatively more massive. 

Most recent analyses of the “autocorrelation lengths” (see Sec. 5) in the rings from 
Cassini data taken at face value seem to indicate that the typical size of particles are smaller 
(implying the opacity is actually higher) in the plateaus than away from them in the low-r 
regions (Colwell et al, 2011; 2012). This could mean, for example, that the largest particles 
in the plateaus may be smaller than the largest particles outside of the plateaus, something 
that would likely affect the viscosity there as well. On the other hand, Marouf et al. (2012) 
find from direct inversion of the scattered Cassini RSS signal that the largest particles are 
much larger in the plateaus than in the surrounding C ring supporting Eq. (22). In this Erst 
paper, we merely explore some implications of our assumed opacity profile with the noted 
caveats until further analyses are conducted that resolve the apparent particle size paradox 
in the plateaus. 


2.2.7 Ring viscosity 

Apart from the angular momentum transport due to the ballistic transport process, the other 
non-Keplerian effect we consider in our model is viscous angular momentum transport. The 
inclusion of viscosity v as a diffusive mechanism is key in counterbalancing the sharpening 
effects of ballistic transport, which can lead to unrealistic growth of (e.g., plateau) edges. Well 
prior to Cassini, the efforts to quantify ring viscosity both analytically (e.g., Goldreich and 
Tremaine, 1978; Araki and Tremaine, 1986) and numerically (e.g., Wisdom and Tremaine, 
1988) assumed ensembles of identical spheres undergoing inelastic collisions while orbiting in 
the gravity held of the planet. A great advantage of these models, despite various simplifying 
assumptions, is that the description of the kinematic viscosity is expressible in a relatively 
simple, analytical form which Durisen and colleagues employed in their original structural 
evolution investigations. 

For the bulk of the simulations presented here, we utilize a radially variable expression 
for the ring kinematic viscosity used by previous workers (Durisen et al, 1992; 1996; CE98) 
which is based on the work of Wisdom and Tremaine (1988), whose calculation of the viscosity 
includes both local and nonlocal contributions to the angular momentum flux. Durisen et 
al. (1992) provide an analytical fit to the numerical results of Wisdom and Tremaine (1988) 
for an ensemble of identical particles of size r 


3v 

0 % 


2.5 r 1 ’ 1 + 


1.7 r 

0.30 + r 2 ’ 


(23) 


which has an accuracy of ± 5 — 10% for 0 < r < 2 when compared with the numerical 
results. The first term on the RHS of Eq. (23) is due to the non-local shear stress, while 
the second term on the RHS is due to the local shear stress (e.g., see Schmidt et al., 2009) 6 . 
A suitable coefficient Qr 2 /3 was determined based on the results of simulations with two 
particle sizes (e.g., Salo, 1991) that showed that the overall spreading rate of particles is 
determined by the largest ones. For the r -3 power law that we use, and our fiducial lower 
and upper bounds rq and r 2 (Sec. 2.2.2), the bulk of the mass and 1/3 of the optical depth 

6 Note that we have used the estimate for the velocity dispersion of the rings c ~ 2r£7 appropriate for very 
flattened systems in the local shear stress term (e.g., Daisaka et al., 2001; Scmidt et al.., 2009) 
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are clue to particles with radii greater than r = 64.6 cm (Durisen et al, 1992). Thus the 
expression for the kinematic viscosity in the BT code is 


" (fl ' T) = 0 ' 51 (ed^)' 


1.55 R & 


3/2 r 


1.5 (t/3Y a + 


( r /3) 


0.30 + (r/3) 2 


cm 2 s 1 . 


(24) 


Durisen et al. (1992) found that a steady-state edge like that of the inner B ring can 
be sustained by a balance between BT, with a strongly prograde ejecta distribution, and 
viscous diffusion (using Eq. [24]). This finding was somewhat remarkable, considering the 
large optical depth contrasts between the low- and high-r regions that characterize the inner 
A and B ring edges. 

Choosing different values for r, which effectively amounts to changing the particle size 
distribution upper and lower bounds, allows us to vary the magnitude of the viscosity in the 
code. We will refer to this variation in magnitude for specific model runs through a scaling 
factor f u = (r/64.6cm) 2 which is unity for our fiducial parameters. Along the lines of the 
discussion in Sec. 2.2.6, we expect that we will be able to more reliably estimate f directly as 
we gain a better understanding of the ring opacity from continued analyses of the available 
Cassini occultation data. 

We note that recent analyses of Cassini UVIS, VIMS and RSS occultation data have 
highlighted the “clump-and-gap” structure of the rings, presumably due to the presence of 
self-gravity wakes, with the majority of the ring mass possibly hiding in the opaque clumps, 
while the gaps contain a lower r and perhaps smaller particles (e.g., Colwell et al, 2007; 
Nicholson et al, 2008; Robbins et al, 2010). Consequently, a more robust model for the 
viscosity that takes into account the rings’ wake structure should include the effects of 
mutual gravitational interactions (e.g., Hahn, 2008). We consider such a model in Appendix 
C. Although no wake activity has thus far been observed in the inner B ring and C ring 
(though see Colwell et al, 2007; 2009), we find that the C ring plateaus may marginally 
satisfy the condition for wake formation and we consider this briefly in Sec. 3.3. 


2.2.8 Ejecta yields and velocity distributions 

CE98 assumed the mass ejecta yield for primarily silicate material impacting icy particles 
to be in the range Y ~ 10 4 — 10 5 . Given this range for the yield, the results of Lange 
and Ahrens (1987) imply a value of Y — 3 x 10 4 at an impact velocity of 14 km s _1 , 
which is the value that CE98 adopted for their simulations. However, CE98 pointed out 
that their adopted values are scaled from experiments at lower impact velocities and refer 
to macroscopic projectiles hitting solid icy targets, rather than micrometeoroids impacting 
ring particles, where the regolith grains may typically be the same size or not much bigger 
than the impactor. Because ring particle bombardment continues to be a more complicated 
process than modeled in laboratory experiments, we choose to use the yield as a variable 
parameter with a range from ~ 10 4 — 10 6 , which covers hard to soft, fluffy or porous targets. 

Our ejecta speeds are parametrized by the dimensionless variable x and the ejecta velocity 
distribution is given by a pure power law 
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Xb < x < x t ] 
otherwise. 


(25) 


f(x) 


x~ 


where f(x) dx describes the fraction of ejecta mass having speeds in the interval x to x + dx , 
and Xb, x t are the adjustable lower and upper bounds of the ejecta distribution. Durisen et 
al. (1992) and CE98 adopted a value of n = 13/4 which was a fit to the high-end ejecta 
velocities in hard target laboratory data for solid basalt (e.g., see Burns et al, 1984). If icy 
ring particles act more like sand, then a slightly less steep slope may be more appropriate. 
Here, we adopt n — 3 as our fiducial value. We explore different values of n in Sec. 3.1. 

In this paper, we only model a pure power law for f(x) 7 . Because of the improved 
computational ability of our parallelized BT code, we can go beyond Durisen et al. (1992) 
and simulate large radial regions over long periods of time while simultaneously being able 
to resolve structure for very small Xb- Moreover, the BT code can readily be altered to adopt 
velocity distributions which have a more complex behavior. In future work, we will consider 
velocity distributions that go beyond simple power laws and include isotropic and retrograde 
(relative to Keplerian) dominated impact ejecta distributions due to disruptive impacts (see 
discussion in Sec. 5). These may require that we rerun the ejecta distribution calculations 
of CD90 with altogether different sets of conditions. 


3 Ballistic Transport Simulations 

3.1 Comparison to Previous Work 

In this section, we present some tests of the numerical code. A natural starting point for 
demonstrating the new code’s utility is to reproduce results in earlier papers. Durisen et al. 
(1992) modeled the relatively sharp B ring edge and its transition into a low constant optical 
depth C ring and demonstrated that the B ring inner edge (and presumably also the similar 
inner A ring edge) can be maintained at its presently observed width through a balance 
between viscosity and ballistic transport with a prograde-biased ejecta velocity distribution. 
Durisen and colleagues explored a range of impact yields and ejecta speed distributions f(x), 
with different choices for the lower bound Xb, to find a range over which the B ring edge 
could be maintained over long time scales. Here, we focus on a few specific cases choosing 
the same parameter values as Durisen et al. (1992). 

Durisen et al. (1992) scaled their models using a relationship between the surface mass 
density and optical depth at R = 1.52 Rg such that a = a{r = 1 )r where cr(r = 1) = 96 
g cm -2 . This scaling then represents a constant opacity model for the rings with an effective 
value k — t j (T = 0.104 cm 2 g^ 1 , and, along with the assumption that the ring particles 
everywhere can be described by an r~ 3 power law with our assumed rq and rq , is consistent 

7 Durisen et al. (1992) looked at pure power laws as we utilize here, but also looked at a power law with a 
“knee” which was considered more applicable to model the hard target, higher-end speed distribution. The 
main difference between the two distributions that Durisen and colleagues found was that the pure power law 
required somewhat higher yields to dominate viscous transport and produced a steeper steady-state edge. 
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with mostly icy ring particles with density 1.12 g cm 2 . The initial condition for the B ring 
edge models has the optical depth distribution 

0.950 

r(R, 0) = 0.050 + e ( L52 R s -r)/d i ’ (26) 

where the parameter D = 0.0015 Rg gives an edge width of ~ 90 km, comparable to what 
is observed. 

In Figure 2, we show a B ring inner edge model with Y — 3 x 10 5 , a power-law ejecta 
speed distribution with a lower bound Xb = 2x 10” 4 , and a kinematic viscosity with weighting 
factor f u = 1, matching parameters in Fig. 10a, Durisen et al. (1992). For these choices 
of parameters, the gross erosion time is around 4.5 x 10 4 years. The simulation uses a time 
step of At = 0.0025 t^, and the entire run, which corresponds to a physical timescale of 
4.5 x 10 6 years, took 35 minutes of parallel CPU clock time. The black curves correspond 
to selected times during the evolution of the ring structure, while red curves of the same 
style are the corresponding fractional mass of pollutant. Over the first ten gross erosion 
times or so, the ring edge steepens to a roughly steady-state slope. A ramp at the base of 
the edge starts to develop between 9 and 35 tc and becomes quite prominent by the end 
of the run at 99 t^. The ramp’s growth is produced by a small residual between direct 
mass exchanges and “indirect” radial drifts due to ballistic transport in this region, which 
is essentially independent of the ring viscosity (he., |u]j all | 3> |u™ c |) there (Durisen et al. , 
1992). The undulatory structure (Durisen, 1995) that appears in the B ring for higher yields 
(see Fig. 4) is suppressed by viscosity for this value of Y. Also observed further inwards and 
shown in more detail in Figure 3 are undulations that emanate from the base of the ramp. 
These undulations have long wavelengths and may be being forced by the growing ramp 
structure which itself is quite broad. Alternatively, these undulations may come as a result 
of the linear BTI which favors regions of intermediate r (Durisen 1995; Latter et al, 2012; 
2014a, b). As the ramp grows in width, so too does the wavelength of the low optical depth 
region undulations. These undulations are not so easily suppressed because the viscosity is 
much weaker here than in the B ring (c/. Eq. [24]). 

Also plotted in Figures 2 and 3 (red curves) are the fractional pollutant mass in the rings 
for the selected evolutionary times (see Sec. 2.2.3). Initially, the ring begins with no extrinsic 
pollutants anywhere (/ e (i?, 0) = 0). At earlier times, the fractional mass differences between 
the high optical depth B ring and low optical depth C ring are quite pronounced, with a 
sharp transition across the B ring inner edge. As the simulation evolves, the differences in 
composition between the regions retain a nearly constant ratio (roughly between 1.49 Rg and 
1.54 R,,s), but the sharpness of the transition smears out over time due to advective effects, 
an effect found by CE98 in their compositional evolutions. Less polluted material from the 
higher mass density B ring drifts inwards, displacing more polluted material inwards of the 
boundary, as seen in observations (Cooke et al, 1991; Estrada and Cuzzi, 1996). 

The compositional contrast between low- and high-r regions depends on the opacity 
contrast between the regions. The CE98 profile for the opacity given in Eq. (21) would 
produce higher contrast than shown in Figures 2 and 3, but its simple form cannot readily 
be applied to the structurally evolving models because it relies on a fixed position for the 
inner B ring edge. In a structurally evolving model, one would have to redefine Ri and i? 2 
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as the edge evolved, which may occur either by drifting radially inward, and/or by smearing 
due to viscosity. This problem is not encountered in the opacity model given by Eq. (22) 
because k depends on the value of r, not on R. In Figure 3 the most notable thing is that the 
relative radial position of the decrease in fractional mass begins further radially inward in the 
ring with time. This is probably due to two effects: the rings are brightening further inwards 
as time goes on as brighter material drifts in from the edge and darkening is occurring less 
quickly as the ramp grows in mass (and radial extent) relative to regions inside of it. 

In Figure 4, we reproduce the results of a compilation of simulations over 104 to for 
the B ring edge model for various combinations of Y and xj, by Durisen et al. (1992, their 
Fig. 15). The dashed curve corresponds to the case for Y — 3 x 10 5 presented previously 
(for 99 t(j) with the run extended for Eve more gross erosion times. The dotted and solid 
curves represent cases with yields of Y — lx 10 6 and 2 x 10 6 , respectively, and a lower bound 
in the ejecta velocity distribution of Xb = 1 x 1CF 4 which corresponds to a velocity of ~ 2 
m s _1 . These runs also have the number of bins increased to 1500 from the 1000 used in the 

Y = 3 x 10 5 case, as was done in Durisen et al. (1992). Furthermore, the width parameter 
(see Eq. [26]) is chosen to be much narrower, D = 0.00075 in order to account for the finer 
grid being used. Also plotted in Fig. 4 is the Voyager PPS optical depth profile (opening 
angle B ~ 28°) for comparison 8 . The main difference between the two latter cases is that the 

Y — lx 10 6 case achieves a steady-state sharpness with an absence of high-r undulations, 
while increasing the yield by only a factor of two quite readily produces nonlinear growth 
of undulations in the inner B ring. The corresponding red curves are the fractional masses 
for each run. The highest yield case has the least amount of darkening which may seem 
counter-intuitive, but this is because the actual exposure time is much shorter for the same 
value of t,Q. For the Y = 10 6 and Y = 2 x 10 6 runs, the physical time scale for exposure is 
1.4 x 10 6 and 6.8 x 10 5 years, respectively. 

From the standpoint of structural evolution, we find no significant differences between 
our simulations and those of Durisen et al. (1992). The differences that do appear upon close 
inspection may be clue mostly to order x errors which would be most prominent near the 
sharp edge. We note that the simulations of Durisen et al. (1989; 1992) used an erroneous 
expression for the angular momentum of an ejectum h (Eq. [22], Durisen et al, 1989) that 
propagated to the expression for the function A{x, 6 , 0) (see Sec. 2.2.5) used in their code 
(Durisen et al, 1996). Here we use the corrected expression for A that is derived using the 
appropriate h. Furthermore, we ran the higher yield cases at a smaller time step than used 
by Durisen et al. (At = 0.0025 t^). Our runs were done with a time step of At = 0.0005 
and 0.0002 t^ for the Y = 10 6 and Y = 2 x 10 6 cases, respectively. As an example of code 
execution time, the Y — 2 x 10 6 case took ~ 300 — 400 minutes of dedicated parallel CPU 
clock time to complete. 

In Figure 5, we ran models with the same parameter choices as in Fig. 4, except that, 
instead of constant opacity, we use the variable opacity from Eq. (22), in which the particles 
are assumed to be larger in the higher-r regions. The magnitude of the opacity in the C ring 
and B ring is comparable to what was used in CE98 (specifically, k(t = 0.05) = 0.05 cm 2 

8 The Voyager PPS optical depth was used by Durisen and colleagues, and CE98 in their work. In this 
paper, we also use the higher-resolution Cassini UVIS a— Arae occupation data as our reference r when 
modeling real structure. 
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g -1 , k(t = 1) = 0.009 cm 2 g _1 ). The higher C ring opacity translates to more mass in the 
plateaus compared to constant k, but about a factor of five less mass overall in the C ring 
than in the constant opacity model, because there is much less mass outside of the plateaus. 
The inner B ring mass remains comparable. 

The greater contrast in mass density between in the inner B ring and the non-plateau 
regions of the C ring leads to two clear differences. First, the difference in darkening between 
the C ring and B ring is much more pronounced, as one would expect for the larger mass 
density contrast. In other words, the C ring darkens much more quickly relative to the B ring 
for a given exposure time. The greater contrast in mass also affects the growth of undulations 
in the inner B ring. Whereas in Fig. 4 we see that the Y — 2 x 10 6 case produces large 
undulations, which continue to grow to large amplitude at longer times, the introduction 
of the variable opacity model mutes the growth of these undulations. One reason for this, 
at least in part, has to do with the fact that the lower surface mass density in the C ring 
decreases the amount of material BT throws into the B ring. However, longer term runs 
not shown here do not lead to further growth of the undulations as they did in the constant 
opacity case; instead they remain roughly constant in amplitude with time, somewhat like 
what is seen in Fig. 8. Thus, a second reason may be that the variable opacity itself acts as 
a damping mechanism to regulate growth. 

This effect can be seen more clearly by comparing dr/dt for the two opacity models 
as shown in Figures 6 and 7 (cf. Durisen et al., 1992; 1996). In these figures, the top 
panel shows the individual contributions to the change in r with time due to the direct mass 
exchanges between neighboring ring regions (dotted curve), and due to the BT (indirect) 
advective term. The bottom panel shows a comparison between the combined BT effects, 
the viscous effects, and their sum shown on an expanded scale. The net change (all effects 
included) is given by the red curve. Outwards of 1.52 Rg, which marks the inner B ring edge, 
the BT effects are weak at earlier times and generally are balanced out by viscosity. As time 
progresses and the effects of exposure begin to accumulate, the competition between BT and 
viscosity in the inner B ring become more and more imbalanced, and a net growth begins to 
occur leading to the undulations we see in Fig. 4. It should be noted that the undulations 
seen in the B ring, much like the in C ring, could be arising due to the BTI with the growth 
rates being small because the instability is near marginality (see Latter et al ., 2012). 

Even though the viscosity is strongest inside the B ring edge, the perturbed growth is 
more or less “stationary”, which allows for the undulations to grow “in place”. This is in 
contrast to the C ring where, although the magnitude of the variation in dr/dt is greater 
than what is seen in the inner B ring, the variations are transient (i.e., not stationary) in 
nature and tend to somewhat average out. This is consistent with recent work by Latter et 
al. (2014a, b) who show that B ring wa vet, ra in s have very small wave speeds, whereas the 
C ring waves propogate much faster. Finally, in this region, the viscosity has little effect so 
the solid red and black curves overlap. On the other hand, and as also seen in Fig. 5, the 
variable opacity model in Fig. 7 is characterized by much smaller net effects away from the 
sharp edge. 
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3.2 Extending Previous Work and Global Simulations 

3.2.1 Long Term Stability of the B ring Inner Edge 

In Figure 8, we have simulated the inner B ring edge with constant opacity using a yield 
of Y = 1.5 x 10 6 with Xb = 10 -4 , and f u = 1 for a timescale of 800 to This run is an 
extension of Durisen et al. (1992, their Fig. 15) with the exception that ring torques and 
mass loading are included. With these parameter choices, the absolute time scale of the run 
is ~ 7 x 10 6 years. We have plotted several curves in 100 to intervals to demonstrate quite 
clearly how BT is able to maintain the sharpness of the inner B ring edge over long periods 
of time, a result that was implied by Durisen and colleagues, but not shown explicitly in 
their work. The development of this steady-state edge due to a balance between viscous 
spreading and the sharpening effects of BT for the initially narrow edge and prograde ejecta 
distribution we use is dominated by the direct mass effect (Durisen et al, 1992), as shown 
in Fig. 6. Viscosity, however, does not play a role in the formation of the ramp which is due 
to advection overbalancing the direct mass term (see Appendix B of Durisen et al, 1992). 

Although the overall slope of the ramp region in this simulation appears to be consistent 
with the linear ramp seen in the observations (dash-dotted curve), a “hump” structure forms 
there between 1.51 — 1.52 R5. At this time, it is not clear what the cause of this feature is. 
A possible explanation is that the hump may be a manisfestation of the BTI, but operating 
against a very nonuniform background where it is less well understood. We consider this as a 
topic worthy of further study. In the inner B ring, the undulations appear to be approaching 
a steady-state away from the edge for this choice of parameters. Recall that there are none 
for Y = 10 6 and that they grow to large amplitude for Y — 2 x 10 6 . This result highlights 
the sensitivity of the simulations to the choice of parameters for the constant opacity case, in 
particular. As we have demonstrated previously, the variable opacity model is less sensitive 
to the parameter choices and thus may allow for a more robust parameter space to work 
with, although we again emphasize that the k(t) we use is heuristic with significant caveats 
(Sec. 2.2.6). We have done the complementary case with a variable opacity for the same 
set of parameters. We show only the profile after 800 tc (red curve, Fig. 8). The main 
difference is that the inner B ring edge slowly spreads, which implies that a slightly larger 
yield will be needed to maintain it. We also find that even though the C ring has five times 
lower cr overall than the constant n case, the hump structure in the ramp is slightly lower 
in magnitude and the wavelength of the C ring undulations are longer. The C ring optical 
depth is also larger due to the spreading of the B ring inner edge. Lastly, the undulations 
in the inner B ring are muted, but this may also be due to the lack of a steady-state edge in 
addition to fewer ejecta being delivered from the C ring to the inner B ring. 

3.2.2 Steepness of the Ejecta Velocity Distribution 

As shown in Figure 9 for a yield of Y = 3 x 10 5 , Xb = ICR 4 and /„ = 0.125, we have run 
a set of simulations varying the steepness n of the ejecta velocity distribution. The inner B 
ring is modeled with a larger r ~ 1.5, consistent with the o-Arae UVIS optical depth profile, 
dotted curve, because the simulation for n — 2 (blue curve) leads to the inner edge becoming 
quite sharp and the ensuing undulations in the inner B ring grow to large amplitude for the 
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case where r ~ 1 in the B ring. A model for the fiducial choice of n = 3 is given by the black 
curve. There are few things of note that we can identify in these examples. For higher n a 
“notch” appears in the inner B ring edge somewhat similar to what is seen in the data. We 
also see that steeper velocity distributions associated with larger n values lead to a gradual 
“spilling over” of material from the inner B ring edge, even as the slope of the edge closer 
to its base remains similar to its initial steepness. It is possible to interpret the occultation 
data as showing multiple notches in the edge, with a more gradual drop off in r (dotted 
curve) similar to what is seen for the n = 3 — 5 curves. The spreading of the edge occurs 
because the steeper velocity distributions concentrate more material at the smallest x thus 
weakening the larger-scale effects of BT relative to viscosity, particularly at the higher r. 

We suggest that the notches seen in these simulations are related to the BTI, but un- 
der very nonuniform conditions. The notches appear to be more pronounced for steeper 
speed distributions. The reason for this is that higher n values places more ejecta in lower 
re— values which reduces the ejecta contribution to ramp production. Because the ramp is less 
prominent for higher n, the ability for BT to overwhelm viscosity at the smallest x — scales is 
increased (similar to increasing the // parameter, Latter et al, 2012; 2014a, b) which allows 
for BTI-like behavior. Figure 9 further suggests that the BTI-like growth of waves could be 
producing the hump structure seen in Fig. 8. The hump is at an intermediate length scale 
because viscosity is able to control the smaller scale BTI-like waves for n = 3. 

What may be more interesting is that lower values of n lead to a larger and better formed 
ramp, because as n decreases, more ejecta is shifted to intermediate and higher values of x 
so that BT cannot compete well at small scales. This is also the reason why the inner edge 
sharpens further. We note that in Eq. (15) of Durisen et al. (1992), they also predicted a 
linear ramp for a value of n = 2. In their analytical model for ramp production, the first 
term of Eq. (15) predicts what is seen in Fig. 9 quite well, but the second term becomes 
divergent very close to the inner edge, which we do not see as-of-yet over the simulation 
time of 100 tc- hi any case, the real ring edge does not look like the n — 2 profile, but 
more like a hybrid of low— and high — n results. That is, the real inner edge seems to exhibit 
features characteristic of multiple ejecta distributions with different n — values at different 
r’s and x— values, and with prograde and retrograde symmetries. This strongly hints at the 
considerable complexity involved in the ejecta distribution function which requires further, 
in-depth study. 

3.2.3 Hidden Mass in the B ring 

The surface mass density in the opaque regions of the B ring is not known, but it has been 
pointed out that, if the B ring is an order of magnitude more massive than the current 
estimate, then it could more easily resist becoming polluted and darkened, and perhaps be 
as old as the Solar System (Esposito, 2008; Charnoz et al, 2009; Robbins et al . , 2010). 
To explore this possibility, we show a low-resolution simulation of two hypothetical B ring 
structures over 100 tc; using a yield of Y — 3 x 10 5 , Xb = 10~ 4 and variable opacity in Figure 
10. We have chosen to use the kinematic viscosity, but modified it in such a way that the 
viscosity is able to maintain both the inner and outer B ring edges. The reason for this is 
that straight application of the wake viscosity (Eq. [C-l], Appendix C) to the B ring requires 
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unrealistic yields to prevent spreading of the inner edge and smearing out of structure even 
in the high — a regions (see Sec. 2.2.7). In order to maintain the sharpness of both the inner 
and outer B ring edges, we have incorporated a simple linear increase in v(R) outwards 
of 1.6 R,g , up to a factor of two larger v than what it would otherwise have been at the 
outer B ring edge. The dashed black and red curves are the structure and fractional mass of 
pollutant, respectively, for a “lower bound” of the B ring mass estimated by Robbins et al. 
(2010, u ~ 400 g cm -2 in the core) which is already larger than traditional “Mimas” mass 
estimates indicated by the solid curves. The corresponding “high mass” model (u ~ 1000 
g cm -2 in the core), as suggested by Esposito (2008), is shown by the dotted curves. These 
two cases produce a measurable difference in the amount of darkening that would occur in 
the B ring core, both between the traditional mass model (solid curves), and with each other, 
as well as qualitative differences in composition between the “troughs” in fractional mass 
compared with the rest of the B ring. Likewise the difference between the B ring core and 
the C ring (and Cassini division) emphasizes the notion that if the B ring is indeed more 
massive, then the low optical depth regions would have become overpolluted over the age 
of the Solar System, unless they are somehow either regenerated and/or replenished with 
bright material. In future work, we will be addressing these issues in further detail. 

3.3 Simulations of Observed Structure in the C ring 

The simulations presented in Section 3.1, and those of Durisen and colleagues, used as initial 
conditions simple local models for ring initial structure, namely a moderately sharp edge 
bracketed by low and high constant optical depths emphasizing inner edges such as those 
of the B and A ring. Simulations such as these are useful for determining how some of the 
structure we see at or near these inner edges may be produced by BT. Here we explore 
simulations in which the initial conditions consist of the observed optical depth profiles for 
parts of the C ring in an effort to understand what role BT may play in maintaining observed 
structures such as plateaus over long periods of time. The Cassini UVIS a-Arae occultation 
(Colwell, priv comm) optical depth profile we use as our initial r has been modestly smoothed 
using a three-point boxcar routine. All simulations in this section utilize the pure power-law 
velocity distribution with n — 3 described in Sec. 2.2.8 with ay = 1 x 10~ 4 and ay = 4 x 10” 3 , 
the kinematic viscosity (Sec. 2.2.7) with f u = 1, and our variable opacity model (Sec. 2.2.6), 
unless otherwise noted. 

In Figure 11 we show a suite of simulations for yields of Y = 10 4 (solid curves), 10 5 
(dotted curves) and 10 6 (dashed curves) for a period of 2 to over the bulk of the C ring 
and inner B ring. The radial resolution is ~ 6.5 km per bin. The gross erosion time scales 
associated with these yields derived from Eq. (5) are tc — 1.4 x 10 4 — 10 6 years with the 
longest time scale corresponding to Y = 10 4 , but as mentioned in Sec. 2.2.1 this timescale 
is normalized to r = 1. A representative range of r for the plateaus in the C ring from the 
optical depth profile we utilize (blue curve) is r ~ 0.5 — 0.25. Using our variable opacity 
model, these values then imply that the gross erosion times are ~ 2 — 6 times shorter than 
the values given above for the lower r regions. 

In Fig. 11, the black curves correspond to the structural evolution, while the red curves 
correspond to the composition. It is evident that, for the parameter set chosen, the effects 
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due to BT for a relatively small yield of 10 4 are insufficient to maintain any structure against 
spreading by our fiducial kinematic viscosity. Inside ~ 1.5 Rs, the outer edges of the plateaus, 
which should be maintained if the prograde ejecta yield were large enough, smear out as 
quickly as the inner edges. Also, outside of this location, the integrity of the B ring inner 
edge is lost and smears almost as vigorously as the lower optical depth plateaus. This 
behavior also seems to characterize a yield of Y = 10 5 , although less so. Given that the 
outer edges of the plateaus are becoming less sharp, it seems that a longer term simulation 
would lead to smearing out of most of the structure here as well. On the other hand, a yield 
of 10 6 actually leads to a BT influence that is too large for the plateau regions, because the 
prograde ejecta distribution sharpens the plateau outer edges more than is observed. The 
observed plateaus outside of R = 1.46 R5 have a much more gradual outward increase in 
r (similar in slope to the ramp itself). If Y is really this high, stronger advective and/or 
viscosity effects may be needed to prevent as much pileup of material at outer edges as 
the simulations show. However, the great sensitivity of the effect to Y suggests that other 
properties may come into play. For instance, maybe there is a feedback on Y itself, as the 
particle properties change clue to structural evolution. 

Also shown in Fig. 11 are the corresponding compositional profiles for 2 tc- We note that 
all three compositional evolutions have been scaled to the Y = 10 4 values for the fractional 
pollutant mass f e for the purpose of comparison. The calculated amount of pollutant in 
the higher yield simulations is 10 or 100 times less for Y = 10 5 or 10 6 , respectively, by the 
definition of (Eq. 5), because the amount of extrinsic material that has hit the rings is 
also 10 and 100 times less in these cases. Much like the structural evolution, the evolution of 
fractional mass of pollutant for a yield of 10 4 leads to a profile of composition that is quite 
diffuse, and the features drift inwards as they smear. The compositional profiles for Y = 10 5 
and 10 6 look remarkably similar even though their structural profiles do not. This suggests 
that compositional gradients associated with higher-optical depth structure embedded in a 
lower optical depth background have a more gradual radial variation than those in optical 
depth and, moreover, are insensitive to fine-scale structural details. In all three cases, the 
compositional transition between the inner B ring to the C ring is broad, even though the 
transition in brightness and optical depth between these ring regions is abrupt (Estrada 
and Cuzzi, 1996; Estrada et al, 2003; Cuzzi et a/., 2002; 2009). This is consistent with 
the earlier finding by CE98, using fixed structure, that a broad observed color transition 
across an abrupt optical depth transition was produced by ejecta exchanges. The results 
of Fig. 11 can be thought of as a generalization of this B-C boundary effect in that, at all 
abrupt plateau edges, the pollution profile is smoother than t(R ) and insensitive to fine- 
scale structure. CE98 ascribed the B-C compositional gradient to inward advection of less 
polluted B ring material. 

These simulations suggest that a proper balance between the sharpening effects of BT 
due to a prograde ejecta distribution versus the broadening effects of viscosity could possibly 
maintain outer plateau edges at their current level of steepness over long periods of time, 
but in a rather parameter-sensitive way. The fact that the inner edges of the plateaus do 
not sharpen as strongly in the current simulations argues for a significant component of the 
ejecta distribution to be retrograde relative to Keplerian, perhaps due to disruptive impacts. 
Simply from symmetry arguments, this would have the effect of sharpening the inner edges 
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because a significant component of radial drift could be outwards (from inside the plateau). 
A retrograde component to the ejecta distribution may also play a role in ramp formation 
(Durisen et al, 1992), potentially leading to suppression of the hump that we see forming in 
the ramp in our simulations. 

The observed optical depth profile for the rings (Fig. 1) may indicate that retrograde 
ejecta play a significant role in the C ring and are even consistent with this role increasing 
inwards (Sec. 1). In future work, we will formally include different ejecta distributions 
(see Sec. 5 for a discussion of retrograde ejecta). Durisen and colleagues also found high 
yields to be necessary to maintain the inner B ring inner edge, even assuming only kinematic 
viscosity (see also Fig. 8). The explanation for how BT sharpens an inner edge, where a 
large optical depth contrast exists between the low-r and high-r regions, is complicated, 
even in linear analyses that simplify the process, yet still preserve the essential features (see 
Durisen et al, 1992, for a full discussion). The important thing to note here is that the inner 
B ring edge, and by analogy the inner A ring edge, both satisfy the conditions for optical 
depth contrast that lead to maintaining a sharp, steady-state inner edge, whereas the C ring 
plateau inner edge optical depth contrasts do not. Thus having a retrograde component 
to the ejecta distribution may be essential for maintaining the C ring structure we observe. 
This component, if due to disruptive impacts, may have a very different velocity and angular 
distribution from the cratering impact ejecta modeled in this paper. 

In Figure 12 we show an example of the effects of kinematic viscosity in maintaining ring 
structure for different choices of the ejecta yield. The justification for this parameter study is 
that the bulk of the optical depth may be due to particles of size f, which are different from 
our fiducial choice (f u = 1, Sec. 2.2.7), as might come about, e.g., if the size distribution 
in the rings were narrower or broader than we assume (Colwell et al, 2009; Marouf et al, 
2012). This may be particularly true within the plateaus (see discussion, Sec. 5). Figure 
12 shows three different simulations of the plateau features centered at R = 1.465 Rs using 
yields of 10 4 (black curve), 10 5 (red curve) and 10 6 (blue curve) over a period of 2 t a for 
three different choices of f u such that Y / f v — 10 5 (see below). This means the red curve 
corresponds to the Y = 10 5 curves in Fig. 11. The initial optical depth profile is given by 
the dotted curve. The model runs with f u — 0.1 and f u = 10 translate to r = 20.4 and 204 
cm, respectively. The characteristic optical depth for these plateau features is r ~ 0.25 — 0.3 
so that the gross erosion times are about ~ 5 times shorter than the timescales given for 
Fig. 11. 

One point that the parameter choices in these simulations are meant to demonstrate is 
that their exists a degree of scaling between f u and Y which can produce the same set of 
results through variation of these parameters (c/. Latter et al, 2012, their [i parameter). 
This was pointed out by Durisen and colleagues who conflated Y and v into an “effective 
yield” in their calculations. As long as the functional forms of t(R, 0), u{R,t), f(x), and 
<t(t) remain fixed, each simulation represents an infinite family of solutions (Durisen et al, 
1992). The reason for this is that, if these functional forms are fixed, the changes with time 
depend only on a dimensionless ratio of the ballistic transport time scale to, to the viscous 
timescale t u = L 2 / 3z/ where L is the length scale over which viscosity acts. With the variable 
opacity model, however, it is not clear that this scaling invariance will hold over very long 
time scales. Figure 12 indicates that it is still valid over short time scales. 
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Thus in all three cases, even though the absolute time scales between them range by two 
orders of magnitude between Y = 10 4 — 10 6 , the profiles look similar with the outer edge 
slopes being comparable in steepness, while the slope within the outer part of the plateau 
itself remains similar to the initial one. The fiducial value of /„ — 1 and a yield of 10 5 seems 
to be adequate in maintaining the sharpness of the outer edge, although some degradation 
of the edge, and more notably the decrease in the optical depth of the outermost plateau 
at ~ 1.486 Rg, would suggest at slightly higher value for Y. However, for a yield of 10 4 
and the fiducial /„, the BT mechanism is too weak to maintain structure against viscous 
spreading, and all structural features are smoothed out in only 2 tc- On the other hand, 
a yield of 10 6 and f v — 1 leads to the largest pileup of material at the outer edge as BT 
dominates the viscosity producing a spike in the amplitude of r there (see both the Y = 10 4 
and 10 6 cases in Fig. 11). Although not shown here, for this high a yield structure begins 
to form within the plateau itself; yet, this same combination of parameters is effective at 
maintaining the sharpness of the (larger optical depth contrast) inner B ring edge. This may 
imply that particle size distributions and/or particle properties (which may affect yields) 
are different for the C ring plateaus compared to the inner B ring. The similarity of these 
simulations suggests that the wide range in possible impact yields could be accomodated for 
if the dominant particle size can be different from our assumed value (r = 64.6 cm). More 
work will be required to determine particle properties in the plateaus before we can draw 
any definitive conclusions. 

The main difference between these simulations is in how much the inner edges have 
diffused, with radial diffusion being more evident in the lower yield case where the “real” 
time is longer (for Y = 10 4 , 2 t q = 2.7 x 10 6 years). If we were to run these simulations (as 
well as those in Fig. 11) for longer times t > 2 to, the inner edges would smear out for any 
value of f u because there is currently no mechanism in the code to sharpen inner edges of 
low-to-moderate optical depth. The ejecta distribution we employ at present is prograde so 
that BT can only sharpen outer edges for these low optical depth features (see Sec. 5). 

In Figure 13, we vary the lower bound X), of the ejecta velocity distribution to illustrate 
the sensitivity of the rings’ structural and compositional evolution to this parameter for 
Y = 10 5 and a value of f u = 0.3 which correpsonds to a dominant particle size f = 35.4 
cm. We chose x t = 4 x ICR 3 such that the maximum ejecta velocity n e j = 100 m s^ 1 at 
1.5 Rs. The structural evolution over 2 R; is plotted in the upper panel for ay = 10" 4 
(black curve), Xb = 5 x 10 -5 (blue curve), 2 x HR 5 (red curve) and Xb = 2 x 10~ 4 (green 
curve) which correspond to lower bound velocities of 2, 1, 0.4, and 4 m s _1 , respectively. In 
some ways, lowering Xb seems to have a similar “softening” effect as increasing the viscosity, 
in that the amplitude of the plateau outer edges and the amount of structure within the 
plateaus themselves decrease with smaller Xb ■ However, while the amplitude does decrease, 
the sharpness of the plateau outermost edges only decreases slightly with smaller Xb for this 
choice of the viscosity, although the location of the edge is systematically different for the 
three cases. This result is straightforward to interpret. For an ejecta velocity distribution 
with a slope of n = 3, most of the mass is ejected at Xb and falls off steeply at larger x. The 
throw distance of this low-speed ejecta is small, meaning that the angular momentum that 
it carries is being deposited at locations whose angular momenta are not too dissimilar, and 
the smaller mismatch produces smaller advective radial velocities. 
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The compositional evolution associated with this simulation is presented in the lower 
panel. The fractional mass of pollutant f e is plotted with the initial optical depth profile 
(dotted curve) shown for reference. The most interesting thing to note is that the lowest Xb 
case has the most pronounced compositional variation, whereas the opposite is true in the 
structural evolution plot. The interpretation here is that the pollution transport is much 
more localized for lower ejecta velocity. This is illustrated by noting the amount of the 
pollution in the low— r regions between the plateaus is larger because brighter material from 
the high— r regions has not had enough time to advect into the low— r regions. The largest Xb 
(green curve) emphasizes how a large lower bound can lead to significant pollution transport 
even though the sharpening of the outer edge structure is the largest in this case. A more 
localized transport due to the lower ejecta velocities may allow the plateaus to better retain 
their compositional identity relative to surrounding regions for longer periods of time. 

Finally, in Figure 14 we use a combined viscosity model that replaces our usual kinematic 
viscosity (Eq. [24]) by the Daisaka et al. viscosity (Eq. [C-l]) when the condition for wake 
formation (Eq. [C-2]) is satisfied in the plateaus (see Appendix C). Three curves are shown 
for Y = 10 6 that focus on the plateaus just inside the ramp centered at 1.495 R 5 . The 
black curve corresponds to the simulation from Fig. 11 for the same yield. We have already 
indicated that such a large yield for the plateaus leads to BT that overpowers the kinematic 
viscosity for the nominal f v = 1 , and stability requires a larger particle size f to dominate 
the optical depth (cf. Fig. 12, f v = 10). On the other hand, the red and blue curves, which 
correspond to ring particle densities p of 0.5 and 0.9 g cm -3 , respectively, suggest that wake 
formation in the plateaus may increase viscosity enough to suppress the tendency of BT to 
produce large amplitude ‘spiky” structures at the plateau edges for such large Y. The main 
problem that exists with using the Daisaka et al. formulation at face value is that the B 
ring inner edge does not remain sharp and rapidly diffuses even over this relatively short 
time scale in simulations (not shown) done using Eq. (C-l). In order to overcome such a 
strong viscosity in the inner B ring requires possibly unrealistic yields (Y >10'). If wakes 
are present, the viscosity may indeed be different in the plateaus compared to the lower 
optical depth background, but it remains to determine the nature of viscosity in the inner B 
ring. Overall, Fig. 14 (and Fig. 12) may simply suggest that plateaus are better sustained 
if the plateau viscosity is larger than the canonical kinematic viscosity case used for some 
reason in our models, and the inner B ring viscosity is more in line with expectations for the 
kinematic viscosity. Marouf et al. (2012) find a population of 20 m radius particles in the 
plateaus, for instance, which needs to be reconciled with UVIS data that suggests smaller 
local particles (Sec. 2.2.6). 


4 Discussion and Conclusions 

We have introduced a modernized evolutionary code for the modeling of micrometeoroid 
bombardment and ballistic transport in planetary rings that is capable of simulating both 
structural and compositional changes simultaneously over a broad range of spatial and tem- 
poral scales. Our code is based on the original work of Durisen and colleagues (Durisen et 
al, 1989; 1992; 1996), where they modeled the evolution of ring structure with time due 



to BT, and the pollution transport model of CE98 in which structure was held fixed and 
only the relative amounts of non-icy and icy constituents were allowed to evolve due to the 
effects of BT. A major innovation has been the parallelization of the code, which not only 
allows us to explore a broader range of parameter space but will eventually allow for the 
implementation of different ejecta distributions and the ability to treat inhomogeneous and 
non-axisymmetric structure. Furthermore, our code is robust enough to readily incorporate 
new physics (see below). These innovations are important because previous studies on bal- 
listic transport were limited both by computational demands and a lack of key physical data 
to work with. Using new Cassini data, and current computing technology, more detailed 
in-depth investigations are possible. 

Code tests. As a primary test of the code’s utility, we have reproduced some of the 
structural results of Durisen et al. (1992, specifically their Figs. 10a, 10b and 15), who 
simulated a relatively sharp high-r inner B ring edge and its transition into a low, constant 
optical depth C ring using constant opacity. Our results are essentially identical, with the 
exception of some small scale (order ~ x near edges and ~ x 2 away from edges) differences 
in the highest yield cases, traceable to slightly different expressions for the function which 
relates the radius of ejection of an ejectum to that of its possible radii of absorption (Sec. 
3.1). Comparing our parallel simulations to those using a serial version of our code, we found 
that a speedup (using n p = 8 on a quad core CPU) of roughly an order of magnitude. The 
increase in speed over the original serial simulations of Durisen et al. would be significantly 
more than this considering the processor speed of the day. 

All previous simulations assumed the ring opacity n to be constant. We ran the same 
set of simulations as above, but instead used a variable opacity model (Sec. 2.2.6). This 
allowed us to model C ring plateaus which are more massive than would be predicted by a 
constant opacity model. The values for the opacity were scaled to the average values gleaned 
by CE98 for the C ring and inner B ring from occultation-variance results (Showalter and 
Nicholson, 1990). Using CE98 scaling for the opacity in the C and inner B ring leads to a 
surface mass density outside of the plateaus that is about five times less than in the models 
of Durisen et al. (1992) but leads to much more mass in the plateaus (Sec. 3.1, Fig. 5). 
The greater contrast in a between the two ring regions causes considerable differences in 
the amount of darkening that occurs over the same time scale for a given Y, with variable 
K producing a much more polluted C ring relative to the inner B ring. The variable n also 
had a considerable effect on the undulations that are seen in the inner B ring for the highest 
yield case (Y = 2 x 10 6 ), with their growth suppressed relative to the the constant k case. In 
both opacity models, a ramp forms at the base of the inner B ring, although it is not linear 
in its early stages. 

Long evolutions of the B ring edge and formation of the ramp. In order to demonstrate the 
ability of BT to maintain a steady-state inner B ring edge over much longer times than could 
be modeled by Durisen et al. (1992), we have run an 800 t^ case with a yield of Y = 1.5 x 10 6 
(~ 7 x 10 6 years) for a constant opacity. For the chosen set of parameters, we indeed find 
that the inner B ring remains sharp and relatively unchanged over the simulation timescale. 
We also find that the undulations in the inner B ring are approaching a steady state. The 
ramp region continues to develop as brighter material spreads inward over time, with the 
overall slope of the ramp approaching rough agreement with that seen in the observations. 
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However, a hump in the ramp grows as well, in disagreement with observed structure, and its 
presence may contribute to the structure that appears at radii inwards of the hump’s radial 
location. We have run a complementary case using the same parameters for the variable 
opacity. We find that the inner B ring edge is spreading indicating that a slightly larger Y 
is needed in order to maintain it. The main differences, which may be dne to the spreading 
edge, is that the inner B ring undulations are muted, while the C ring undulations appear 
to be longer in wavelength than in the constant opacity case. 

The reason for the formation of the hump is not entirely clear, but we suspect that it is 
a product of the growth of BTl-likc waves at an intermediate lengthscale. Higher values of 
n lead to “notches” rather than humps in the edge that are somewhat similar to the notch 
seen in the observational data (see Fig. 9, dotted curve). This is because a higher value of 
n places more ejecta in lower x — values which reduces the contribution to ramp production. 
Another effect of higher n is for the B ring inner edge to round off more to larger radii, even 
while the slope of the edge near its base remains sharp and similar to its initial value. This 
is because for steeper sloped ejecta distributions, more material is concentrated at smaller x 
which weakens the ability of BT to counter viscosity at high r. For lower n, while the edge 
remains sharp, no hump forms in the ramp region and the ramp appears quite linear (Sec. 
3.2, Fig. 9) in agreement with observations. The wide difference in the effect of the variation 
of n strongly suggests that not just the particle size distribution, but also particle properties 
such as hardness may vary significantly over the rings. In particular, the B ring inner edge 
seems to exhibit features that suggest that multiple ejecta distributions with different r, n— 
and x — values are in play, as well as prograde and retrograde symmetries. 

C ring plateaus. We have run a series of models that focus on simulating how BT 
can maintain the observed plateau structure within Saturn’s outer C ring, studying the 
effects of varying the yield Y , the strength of the kinematic viscosity /„, the possible role 
of gravitational wakes, and the lower bound of the ejecta velocity distribution X),. For our 
fiducial model parameter choices, the sharpening effect of BT for relatively small yields (10 4 , 
but less so for Y = 10 5 ) is insufficient to overcome the effects of viscosity. Outer edges of 
plateaus, which should sharpen due to BT alone for a prograde ejecta distribution, smear 
out as quickly as inner edges. For the highest yield case (Y = 10 6 ) BT effects more than 
overcame viscosity in the C ring plateaus, leading to oversharpening of the plateau outer 
edges; however, yields > 10 6 are required to maintain the sharpness of the inner B ring edge. 
Thus BT effects are sensitive to parameters, which may vary from place to place. The effects 
of BT on composition globally lead to a broad transition in f e across the inner B ring edge 
and C ring, consistent with what is observed and previous models (CE98). ft is most gradual 
for the lowest yield, which represents the longest “real” time evolution. 

We also varied /„ for yields of 10 4 — 10 6 (while holding xj constant) to demonstrate the role 
of the kinematic viscosity in maintaining ring structure. We find that suitable combinations 
of f v and Y can be found between Y = 10 4 — 10 6 , with larger Y requiring larger /„, that 
produce the best fits to the shapes of C ring plateaus. Agreement is not exact, but there 
is strong evidence for a combined dependence (through the parameter Y/f u ), as also found 
by Durisen et al. (1992). We have also noted that the C ring plateaus marginally satisfy 
the condition for wake formation, calling for a different viscosity model. Using the highest 
yield case, where large increases in particle size f are needed for a purely kinematic viscosity 
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model to prevent BT from oversharpening structural features, the wake model for v inhibits 
oversharpening much more effectively without having to resort to very large values of the 
average particle size f. However, applying such a model to the inner B ring edge then 
requires possibly unrealistic ejecta yields (Y > 10 7 ) to maintain a sharp edge against viscous 
spreading. 

The lower bound of the ejecta distribution defines the minimum throw distance of mate- 
rial resulting from an impact. The lower the ejecta velocity, the less will be the difference in 
the specific angular momentum between where the material was ejected and where it lands. 
This leads to smaller radial drift, even if the amount of material being deposited for a given 
Y is the same. This has two effects on the plateau edges. The first is that the location of 
the edge is markedly different for different choices of the lower bound of the ejecta velocity 
distribution (less overall drift for lower Xb) for the same evolution time scale (Sec. 3.3, Fig. 
13). Meanwhile, for all values of Xb, the slope of the outer edge remains roughly the same. 
Second, lower ay has a tendency to soften the secondary structure within the plateaus for 
cases in which v is weak compared to BT, as well as the “spikiness” of the outer edges. 
Though lower ay leads to less dramatic structural changes, we find that the lower ay also 
produces more contrast in composition, because the direct term dominates the advective 
term, and pollution transport is much more localized. 

Sensitivity to particle properties, ft seems clear from our preliminary simulations that a 
better understanding of local particle properties is needed (he., f, ri, ry, ?>, n, ay, ay, and the 
particle’s bulk density) in order to better constrain the viscosity, yield and ejecta distribution. 
If the particle size distribution is different in the plateaus compared to the low— r regions 
outside the plateaus, or in the ramp or inner B ring, then a viscosity that only depends on 
r (or a) and a constant yield is probably insufficient to describe BT over the large spatial 
regions we are modeling. A lower yield, as would be expected for hard targets, can match the 
structure of the plateaus if the viscosity is also systematically lower compared to our fiducial 
choice ( f u = 1). These values may very well need to be significantly different in the inner 
B ring, where Durisen et al. (1992) found that much larger yields were needed to maintain 
the B ring inner edge for a kinematic viscosity. A wake viscosity is even more problematic. 
Finally, not only is the ejecta distribution likely more complicated than we assume here, it 
is also likely variable across the rings. 

Pollution of a more massive B ring: We presented a toy model for the compositional 
evolution within the B ring, addressing the idea that the B ring may be more massive 
than previously thought, which would have important implications for the rings’ age. We 
compared two massive models suggested by other workers (Robbins et al. , 2008; Esposito, 
2008) to a standard “Mimas” mass model in order to demonstrate both quantitative and 
qualitative differences in darkening over time that would occur in the B ring core, as well 
as the overall contrast between the B ring and the C ring. Our simulations (Fig. 10, Sec. 
3.2.3) only span 100 t^ (~ 5 x 10 6 years), but if a ^ remained constant since the LHB and 
we extrapolate the results to cover this length of time, then the amount of darkening in 
the B ring core would be a fraction of a percent for the high end mass model, which is 
consistent with the work of CE98. However, the rest of the rings, especially the C ring and 
Cassini division would be far more polluted than the observations indicate. In all likelihood, 
the micrometeoroid flux was even higher early on, which would only serve to heighten the 
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compositional differences between low and high r regions. 


5 Future Work 

The simulations we present here, along with the work of Latter et al. (2012; 2014a, b) strongly 
suggest that BT (and BTI) is actively at work within the rings. In particular, the formation 
of ramps and the sharpening of edges can be attributed to the effects of micrometeoroid 
impacts and ballistic transport of impact ejecta. Our simulations clearly imply that there 
are still missing physics that will need to be implemented into the code. A challenge for the 
future will be to constrain the microphysics within the rings in order to understand how all 
the various effects of BT can be present at the same time. For example, how can we reconcile 
the form of the C ring plateaus, ramp and the B ring inner edge relative to each other? The 
sensitivity of our results to the variation of k, u, Y and the form of the ejecta distribution 
indicate that simple models in which these properties are considered to be the same across 
the rings are inadequate, and will need to be studied in more detail. Below, we summarize 
improvements that will be applied in future work. 

Isotropic and Retrograde Distributions. Disruptive collisions tend to fragment the tar- 
get ring particle into many pieces, leading to ejecta distributions that differ from the non- 
disruptive (cratering) case in terms of yield, velocities and directionality. The “yield” consists 
of most (or even all) of the target particle, with most of the mass located in the largest frag- 
ments (e.g., Nakamura and Fujiwara 1991) which move at much smaller velocities relative to 
those associated with cratering impacts. Lower relative velocities imply that the throw dis- 
tances of most material in a disruptive impact are smaller, and their ejecta may mostly move 
in the direction of the projectile (retrograde relative to Keplcrian, Paolicchi et al., 1989). 
Thus a considerable fraction of disruptive collisions could result in an ejecta distribution in 
which most of the mass may land at radii inwards of their ejection point. Projectiles which 
only marginally disrupt a ring particle may only result in an isotropic velocity distribution 
(Nakamura and Fujiwara 1991), but even this would be less prograde than the prograde 
distribution we assume for cratering ejecta in this paper. 

Retrograde or isotropic ejecta may play a significant role in evolving local ring structure. 
We have already indicated that difference between the effects of prograde and retrograde 
velocity distributions may be the most important in the C ring plateaus. The characteristic 
outward-increasing slope in r seen in all the plateaus outside 87500 km (Fig. 1) is a well 
known signature of BT due to prograde ejecta. Retrograde (or even isotropic) ejecta would 
lead to the opposite (or flattened) trend, which is seen in all the plateaus inside 87500 km. 
An observation suggesting that disruptive impacts may be playing a role is recent Cassini 
Radio Science (RSS) data that imply there are fewer ~ 1 — 3 cm-sized particles in the 
plateaus than outside the plateaus (Cuzzi et al., 2009). This is significant because ~ 1 cm 
particles may be the size most easily disrupted by meteoroid bombardment, if the peak in 
the micrometeoroid flux mass distribution indeed occurs near a radius of ~ 100 pm. For 
instance, one can crudely estimate the ratio of the mass of the largest fragment to the 
original target mass j\ ~ (e.g., Melosh, 1989) where p ~ 10 6 is the corresponding 

target /impactor mass ratio, and Q* ~ 10 6 ergs g -1 (e.g., Stewart and Leinhardt, 2009) is the 
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fragmentation energy per unit mass. Including gravitational focusing, the impact velocity 
tv ~ 30 km s' 1 in the C ring, giving /*, ~ 0.1, well within the regime of catastrophic 
breakup. Whether or not this suggests that disruptive impacts occur more prominently in 
the plateaus is something we will explore. 

Recent work that uses observations by New Horizons to estimate the flux of EKB dust at 
Saturn suggests that this may be an equally important source of extrinsic material (Poppe et 
al, 2010, Poppe and Horanyi, 2012). Furthermore, though the mass peak for this population 
(a few microns) is likely too small to lead to disruptive impacts, these authors’ models 
indicate that the ejecta distributions produced by EKB dust grain impacts are likely to be 
more isotropic (Poppe, priv. comm.), because their inflow distribution to the rings is more 
isotropic. Moreover, this population is associated with a large gravitational focusing factor 
(CD90). 

Ring opacity and particle properties. The ~ 100 spiral density waves analyzed by Cassini 
UVIS and VIMS (e.g., Colwell et al., 2007; 2009; Hedman et al, 2007; Nicholson et al., 2008) 
have greatly improved our knowledge of key ring dynamical properties such as the optical 
depth and mass surface density, leaving the surface mass density unconstrained only in the 
dense B ring core. Other workers find that, even on smaller scales than the core region, 
the rings are composed of virtually opaque clumps separated by nearly transparent gaps. 
This structure has been attributed to self-gravity wakes, which numerical simulations have 
shown to form spontaneously as a result of the combined effects of self-gravity and collisional 
damping of particles (e.g., Salo 1995; Daisaka and Ida 1999; Ohtsuki and Emori 2000; Lewis 
and Stewart, 2009; Robbins et al, 2010). Nicholson et al. (2008) and Colwell et al. (2008) 
have exploited this expectation to separate the large and unconstrained r within dense wakes 
from the r and fractional area of the nearly transparent gaps. These workers find that, for 
(VIMS) occultations spanning both the A and B rings (11° < B < 51°), gap optical depths 
are relatively constant (0.15 < r < 0.3) wherever wakes exist. 

Analysis of these extensive observations (which also gives the typical pitch angle, height 
and spacing to width ratios of the dense wake clumps, e.g., Colwell et al, 2006; 2007) 
allows one to determine a local autocorrelation length scale which characterizes the physical 
dimensions of contiguous physical material (clumps) blocking the starlight. This scale varies 
depending on the viewing geometry, because the clumps are not spherical, thus complicating 
the search for an “effective largest particle size” that characterizes the ring opacity n. By 
combining results from multiple Cassini UVIS stellar occultations that encompass a broad 
range of viewing geometries, it may be possible to constrain the absolute physical dimensions 
of the self-gravity wakes (Colwell, priv. comm.), which can then be applied to our opacity 
model in the BT code. 

Better understanding of viscosity. According to more recent dynamical simulations 
(Stewart et al, 2007; Robbins et al, 2010), gravity wakes in the A ring behave quite dif- 
ferently than they do in the B ring. While clump lifetimes are much shorter in the A ring 
relative to the B ring, clumps extend further in the radial direction when stretched, and 
the angle of a wake relative to the azimuthal direction varies more with time as it evolves 
(Colwell et al, 2009). These differences will clearly affect the viscosity, but there is still no 
clear indication of how this happens or how the magnitude of the viscosity varies radially 
at different points in the rings. Some W-body simulations indicate the possibility of order 
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of magnitude changes between 120000 km to 100000 km for a constant a (Stewart, priv. 
comm.). In order to apply these updated models with confidence to our BT code, we will 
need a clear idea of the relationship between v, a and r over a broad range of surface density 
values. This is important because the simple model for wakes we have applied in this work 
(Daisaka et al, 2001; Figures C-l and 14) predicts a viscosity in the inner B ring that is up 
to ~ 2 orders of magnitude greater than what is predicted from a purely kinematic viscosity, 
just how much greater depends on the actual r which is a function of elevation angle, and 
also wake orientation if present. In any case, such a large v would make it very difficult 
for BT to maintain the inner B ring edge in steady state. Work in this area is ongoing and 
improvements will be applied to future work as they become available. Finally, other models 
and parameterizations for the viscosity exist (e.g., see Schmidt et ah, 2009) which we intend 
to explore in our models going forward. 

Code, modifications. Our BT code currently uses a rezoning scheme to maintain a rel- 
atively uniform width for the Lagrangian cells over the course of a simulation. If left 
unchecked, bins can become so narrow as to lead to unrealistic structure, or, worse, bin 
edges can even overlap. The process of rezoning as currently done in the BT code avoids 
these problems but requires some averaging, which has the effect of smoothing quantities and 
thus artificially affecting structure. This will likely become much more of an issue for very 
long simulations that may span thousands of t^. The rezoning algorithm will be revisited 
with the idea of developing a method for handling the merging and splitting of bins that does 
not lead to the systematic smoothing of structure. Since we will be interested in modeling 
local (hner-scale) structure as well, we feel this is a prudent endeavor. 

From our results here, it would seem that the ring particle size distribution and other par- 
ticle properties may vary considerably over the rings. Previous work avoided this complica- 
tion because it introduced one or more parameters in the already computationally expensive 
calculation of the losses and gains. However, with the parallelization of the code, the ability 
to expand the integration of the T’s and A’s is no longer a roadblock to further progress. 
In future papers, we expect to include more treatment of variable particle properties as well 
expanding the number of points used in the integration over the velocity distribution x, 
which will also be useful in local, hner-scale simulations, as well as expanding the number 
of points used in our integrations over the azimuthal angle <fi. 

Longer evolutions with improved physical models continue to provide encouragement 
that BT is an important process driving ring structure. The differences in how structure and 
composition behave under BT, and the localized versus more global transport of pollutants 
may be something we can distinguish from the observations. In forthcoming papers, we will 
be modeling both specific features as well as large scale structure based on the observed 
optical depths from Cassini UVIS, while the ratios of Cassini ISS color images can provide a 
measure of the differences in composition across these features. Ballistic transport modeling 
along with a combination of radiative transfer modeling in the spirit of what was done in 
CE98, combined with the very extensive data set provided by Cassini, might then provide 
much tighter constraints on how pollutants are transported within the rings. It ultimately 
may help to answer two long-standing questions: how old the rings really are and what they 
are made of. 
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Appendix A: The ejecta yield function 


We utilize a radiative transfer analog, which assumes plane-parallel symmetry and treats in- 
coming meteoroids and their ejecta as “photons” , to compute the ejecta distribution function 
*3/ which has been parameterized from the numerical results of CD90 in the form: 


SUa.ft R,t) = -T g(R, T )X e ~(°A.,) 8 /s cos (a/2) (cos /3) m . 
ltr ctp cos (ccp/zj 


(A-l) 


The ejecta distribution function is most conveniently defined in terms of the cone and clock 
angles (cc,/3) of the emergent ejecta measured from the ring particle velocity vector and are 
related to the zenith and azimuth angles (0, (j)) of the emergent ejecta through the relations 
cos 6 = since cos f3 and tan (f> = — sin /Stance, respectively (see Figure A-l). The function 
g(R } t ), and the parameters a p , m and S are themselves functions of r, a and/or R (Equations 
[41-46], CD90). The function g(R,t) is used to calculate the probability that incoming 
micrometeoroids, and/or subsequent incoming or outgoing ejecta, will impact a ring particle 
as they pass through the rings. That is, the ejecta distribution 3/ is direction-dependent 
and describes the total rate of impact ejecta emission or absorption per unit area and time 
for a given local ring region of surface mass density a and optical depth r. The distribution 
function above assumes that impacts are cratering and non- disruptive (see Sec. 2.2.8). 
Figure A-2 shows an actual calculation of 3V as a function of a and /3 for optical depths of 
r = 1 and r = 0.1. In these plots, the prograde bias is clearly illustrated by the maxima in 
the distributions occurring for a < 1.5 radians. The parameterized fit in Eq. (A-l) matches 
the actual results quite well (see Fig. 13, CD90). 


Appendix B: Secondary effects to radial drift 

In this appendix, we discuss the secondary contributions to radial drift in more detail which 
are due to direct mass and angular momentum deposition rates. If one considers the mass 
deposition rate alone (“mass loading”, second term on the RHS of Eq. [14]), the timescale 
for a ring to be hit by its equivalent mass is given by t mas i = o-/cXi m = Etc (c/. Sec. 2.2.1). 
A comparable drift time from some radial location R to the central body due to this term 
would be AR/|u“ asl | = (AR/2R)t mas \. For typical ring parameters, radial drifts due to the 
mass loading effect alone would cause the C ring and most of the B ring to be lost to Saturn 
faster than it would take for the accumulation of enough extrinsic material to equal the 
original o (Ip, 1983; 1984; Durisen et al, 1996). 
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The angular momentum deposition rate leads to a net negative torque on the rings, 
and thus an additional radial drift component Up° iq (Ip, 1984). The reasons for this can be 
twofold. The first is that in regions of moderate or low optical depth, micrometeorites can 
pass through the ring without being absorbed. Recall that the reason that the ejecta flux 
distribution of CD90 is strongly biased in the prograde direction is due in great part to the 
aberrations introduced by the motions of Saturn and the orbits of the ring particles. Due to 
these aberrations, the slant path of the incoming micrometeoroids statistically leads to more 
impacts on the leading face of orbiting ring particles, and these particles are more completely 
absorbed. A negative torque occurs because meteoroids with negative specific angular mo- 
mentum are preferentially absorbed, and this would be the case even if the micrometeoroid 
flux were isotropic in the frame moving with Saturn (Durisen et al, 1996). This effect does 
not affect very high optical depth regions because then all micrometeoroids are absorbed. 
A second effect that can lead to a negative torque is that a significant amount of the ejecta 
might escape the system, carrying with them specific angular momentum. With very large 
yields, even a small fraction of ejecta that might vaporize, or obtain large enough ejecta 
velocities could possibly lead to significant changes over very long evolutionary times, but 
we do not consider this possibility currently in our models 9 . 

Durisen et al. (1996) used the techniques of CD90 (e.g., see their Fig. 15) to calculate 
the drift velocities due to mass loading and asymmetric micrometeoroid absorption, which 
causes the negative torque, as functions of optical depth for a broad range of radii and optical 
depths. They derived simple analytical functions of R and r, which we utilize in the BT 
code, that capture the results of their numerical calculations. The analytical fits to the drift 
due to collective effects are given by 

4 orq (i?, r) = - [1.1 x 10 ~ 8 (R/R s ) + 4.0 x 1(T 8 ] e (- T / a28 )°' 74 cm s" 1 , (B-l) 


.masp d 1 .1 X 10 8 (R/Rs) + 6.0 X 10 8 ^ _ e -r/0.47^ 0 ' 98 cm g -l 


u“ asl (R, r) = 


(B-2) 


In these expressions, it is implicitly assumed that the ring mass surface density is 96 g cm 2 
for t — 1 . 


Appendix C: Viscosity in the presence of self-gravity 
wakes 

The ubiquity of gravitational wakes in Saturn’s A and B rings has emphasized the importance 
of the rings’ own self-gravity (e.g., Colwell et al, 2007; Nicholson et al, 2008; Schmidt et al., 
2009; Robbins et al, 2010). The relative abundance of these “clumps and gaps” controls the 
observed optical depth, but the ring viscosity is controlled by the mass and dimensions of the 
clumps and/or large particles. Numerical simulations demonstrate that once strong wake 

9 This does not violate our “thin disk” assumption in Sec. 2.2.2. Even if some ejecta obtain a relatively 
large u e j, the overwhelming bulk of impact ejecta will only achieve velocities in the ~ 1 — 100 m s -1 range. 
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structures form, the effective viscosity is dominated by gravitational torques due to the wake 
structure, because gravitational scattering by wakes increases relative random velocities, 
which enhances the local component of the viscosity (Salo, 1995; Daisaka and Ida, 1999; 
Schmidt et al.,. 2009). 

In the absence of wakes, the kinematic viscosity (Eq. [24]) would dominate. On the 
other hand, introducing self-gravity allows the possibility for transient wake formation, as 
demonstrated in local IV-body simulations (e.g., Salo, 1995; Daisaka and Ida, 1999; Ohtsnki 
and Emori, 2000), if an increase in the local mass surface density sufficiently decreases the 
Toomre parameter such that Q = Qc/3.36Ga < 2. One model which includes this physics 
is that of Daisaka et al. (2001), who derived an analytical expression for the viscosity in the 
presence of wakes that can be readily utilized in our code: 


G 2 rf 0 

u(R,a) = 26-^ 


a 



(C-l) 


where r* h = (. R/2r)(2m/3M ) 1 / 3 is the scaled, mutual Hill radius of ring particles of size r, 
and thus the viscosity depends on the particle internal density. The numerical factor in Eq. 
(C-l) is a fit to their IV-body simulations, and quantities with subscript 0 are evaluated at 
R s . 

In terms of optical depth, the condition for wake formation (Ohtsuki and Emori, 2000) 
can be expressed as 


r > 0.45 



/ 1.45 Rs 

V R 


3/2 


(C-2) 


where po = 0.5 g cm -3 . This expression assumes spherical particles of a single size, and that 
the velocity dispersion for Saturn’s rings is likely determined by the escape velocity from the 
largest particles (Salo, 1995; Daisaka and Ida, 1999; Ohtsuki, 1999). Equation (C-2) suggests 
that the gravitational instability can marginally be satisfied in the C ring plateaus, especially 
if ring particles are not porous. This is better illustrated in Figure C-l where we plot the 
viscosities associated with the two models that we use for the observed optical depth. Here 
we have used the optical depth from the Cassini UVIS u-Arae occultation (elevation angle 
B = 54.4°, Colwell et al, 2009), and assumed a variable opacity (Eq. [22]) with the values 
in the C ring and inner B ring consistent with CE98. In the C ring, the viscosity in the 
low-r regions is dominated by collisions alone, whereas in the highest- r plateaus, as well as 
the inner B ring, the viscosity could be dominated by wakes. 

If the viscosity in the inner B ring is up to two orders of magnitude larger than what 
collisions alone would predict (as we would get using the Daisaka et al. formulation, see 
Fig. C-l), the BT mechanism would require unrealistic yields in order to maintain a sharp 
inner edge. The actual normal optical depth in the B ring outside of wakes depends on 
elevation angle (while in the C ring r is less sensitive, Colwell et al, [2009]). However, wake 
models predict that the observed r has a strong dependence on the orientation of wakes 
relative to the observational geometry and elevation angle, invalidating the classical scaling 
with elevation angle B (Cuzzi et al, 2009) for at least the lower elevation occupations. 
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Nevertheless, even r ~ 1 produces a rather large variation in the wake viscosity relative to 
the highest plateaus due to the a 2 dependence. Thus the circumstances under which the 
viscosity in Eq. (C-l) is applicable is less straightforward to determine (see discussion, Sec. 
5). Furthermore, the possible wake effect is magnified by our expression for k(t) (Eq. [22]), 
which was motivated by the plateaus and may not hold in the B ring. 


Appendix D: Numerical methods 

Stepping in time: In the BT code, the ring is divided into N ‘ringlets’, or annuli with bin 
centers Rj. The annuli are treated as Lagrangian cells whose IV + 1 edges bj ( bj < Rj < b ]+ i) 
drift at velocities due to a combination of the ballistic transport mechanism, ring viscosity, 
mass loading, gain or loss of ejecta momentum, and/or torques that arise as a result of the 
asymmetric absorption of meteoroids (Durisen et al., 1996). The main quantities of interest: 
the surface density er, optical depth r, opacity n, viscosity u, and fractional mass of non-icy, 
extrinsic pollutant f e , are evaluated at bin centers, which lie midway between the moving 
bin edges. Given these quantities at bin center locations R'j at time t n , we advance the 
simulation 


t n+l = f n + Ai n (D-l) 

using a variable timestep A t(t). We describe the execution of our code below. 

1. We implement a rezoning scheme (Durisen et al, 1992) that searches over the radial 
grid to find annuli that have become narrower or wider than pre-assigned limits. If the 
condition is satisfied, the narrowest bin is merged with its smallest neighbor, while the 
largest bin is split in two. The physical quantities of a merged or split bin are adjusted 
accordingly. Typical minimum and maximum limits are 2/5 and 5/2 the mean initial bin 
size, respectively, where the mean bin size is the radial extent of the grid divided by N. 

2. The viscosities u" are updated from the last timestep t n using our chosen viscosity 
law (see Sec. 2.2.7 and Appendix B). 

3. The appropriate timestep size A t n for the current iteration is determined from the 
stability criteria (see below). 

4. The losses A” and gains T n are determined at every if”. Our method of integration 
(described below) uses a lookup table of current values of the fractional mass, the optical 
depth and other quantities that depend on these parameters. Thus, prior to integration, the 
relevant quantities are interpolated onto a very fine grid. 

5. The drift velocity due to ballistic transport at bin center Rj is calculated from Eq. 
(12) using the difference formula 


ball\n 
V R )j 




— A n 
A i,h 



(D-2) 


6. The viscous drift velocity Eq. (13) is calculated at points midway between bin centers 
using the differencing 
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(D-3) 


. . . 3 (avy/R)™ — (auVR)l i 

/^visc'vn v 'J — 1 

) j— 1/2 -n dii Dn pn 

a j-l R j-l R j - R j~ 1 

where all quantities are evaluated at time t n , and the bar represents an average of the quantity 
between spatial bins j and j — 1 (e.g., Rj-i = ( Rj + Rj- 1)/2). The viscous drift velocities 
are then interpolated to the true bin edges which generally do not lie midway between bin 
centers. 

7. If mass loading and ring torques are included in the simulation, then corresponding 
velocities (which depend on r”) are calculated at bin centers using Eqns. (B-l) and (B-2) 
from Appendix B. 

8. The total changes to the ring surface mass density due to direct deposition terms 
(RHS of Eq. [8]), as well as direct changes to the surface mass densities of extrinsic and 
intrinsic materials, are calculated using 


<’*8 = y + At " K<Ui + q, - A"j i 

MS', = (/.)>” + At” + AaJ j , (D-4) 

= [1 - ■(/.)“] + Af [(1 - !)/«,) + Act,] 

where Acr e and Ad; represent the integrals on the RHS of Eqns. (10) and (11), respec- 
tively. The new fractional mass due only to direct terms is then calculated from (/e)”^. = 
(a l n+1 /rr n+1 

we^dir/ u j, dir’ 

9. The “indirect” contributions to changes in a, a e and cr; are due to differential drift. 
Bin-centered velocity contributions to the drift velocity Ur are interpolated to bin edges, and 
then combined to evolve bin edges to their new positions at t n+1 : 


!)” +1 = b] + Af> R )” 


(D-5) 


The new bin centers are then taken to be midway between the new bin edges Rj +1 

(b]ii + !>” +1 )/2. 


10. The divergence or indirect term contribution to the surface mass density on the LHS 
of Eq. (8) then is calculated by modifying the direct term to account for bin compression or 
expansion 


a 


n + 1 


= a 


n + 1 


- 4 ? 


j,dir ^n+1 1 


(D-6) 


where A'j and A'j +1 are the computed areas of the annuli j at times t n and t n+1 . The 
divergence term contributions to cr e and o\ is taken into account by the same technique. 

11. The optical depths r" +1 are calculated using the new values of the ring surface mass 
density and our chosen model for the opacity which can itself depend on the optical depth 
or surface density (see Sec. 2.2.6). 
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Stability and time step restrictions'. The BT code uses a variable timestep A t(t) which is 
determined at each step to satisfy several conditions. The Conrant condition (e.g., see Press 
et al, 1992) that arises from Eq. (8) requires that 


At < min 


M?+i + 



(D-7) 


where the bin width is Wj = bj + \ — bj and an average is used of the total drift velocities Ur, 
which have been defined at bin edges by step 9 above. The condition essentially demands 
that the surface density in the bin being evaluated is not receiving information from beyond 
the boundaries of the spatial regions defined by each grid cell. A second condition that we 
impose is the viscous Courant condition. Code stability requires that the minimum timestep 
be shorter than half the shortest viscous diffusion time across an annular cell, namely 


At < min 


2w* 


|_3(z/j+i + Vj). 


(D-8) 


The viscous contribution to the drift velocity is present in Ur, but the stability criterion 
in Eq. (D-7) is typically less important for simulations with viscosity which we generally 
consider. 

A last condition that the BT code must adhere to is 


At < min 


a' 


(dim )jY 



(D-9) 


which essentially states that the time step must be less than the local gross erosion time of 
the ring. This obviously applies most strongly in situations where the surface density of the 
ring is very low, but not zero, such as near gaps or true edges. Radial cells that are empty 
(<x = 0) are not considered in setting At. At the beginning of each time step, these three 
conditions are checked, and the new At is taken to be the minimum of these, or a default 
maximum time step size which is chosen at the beginning of our simulation. The maximum 
Afmax is 0.00005 — 0.001 t^, depending on the parameter choices for a given run. 

Numerical technique for calculating A ’s and T ’s: Numerical integration of the loss and 
gain integrals (Sec. 2.2.5) represents the most computationally intensive piece of the BT 
code. The integration over 6 is done using a Nq = 21 point Simpson rule evenly spaced 
in cos 9, while integration over f is done using two Gaussian quadrature points over 
the azimuthal angle f = n/2 to tt (the CD90 ejecta angular distribution function possesses 
reflection symmetry through the planes defined by (j) = 0 and f = tt/2, see Fig. A-2, and 
Durisen et al, 1992). Since we primarily use a power- law representation of the velocity 
distribution (Sec. 2.2.8), the ^-integration is done using a trapezoidal rule with typically 
N x — 10 — 20 quadrature points evenly spaced in logo;. 

The value of the 0-integrand at a given Simpson point is determined by interpolation over 
a fine grid of the values for the fractional mass, optical depth, and optical depth-dependent 
functions. Prior to each numerical integration, we pre-interpolate a much finer mesh of points 
Amesh 3> N where the array index of the mesh is chosen to be directly proportional to radial 
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position. This allows for integrands to be evaluated using a mesh value that is closest to the 
actual value of R or R' (Durisen et al, 1989). This procedure somewhat mitigates the costs 
of doing the angular integrals which would require ~ x Ng x N distinct searches and 
interpolations otherwise. Our mesh size ranges from Nm PK h = (50 — 100) x N. Use of the fine 
mesh provides a boost to computational speed at the expense of a small loss in accuracy. 

After other efforts to make the integrations of the losses and gains more optimal and 
efficient, significant improvements in speed are best achieved through parallelization of the T 
and A integration subroutine. This was done using the language independent communication 
protocol Open MPI (Message Passing Interface). Prior to reaching the losses and gains 
integration routine, the ring (composed of N annuli) is divided into n p chunks, where n p is 
the number of processors (“tasks”) being utilized, with each task executing its ‘piece’ of the 
loss and gain integration loops. After the integration is completed, the changes in the values 
of each radial interval must be communicated to all processors before the global calculation 
can continue. That is, outside of the losses and gains routine, the code runs n p separate, but 
identical versions of the simulation until the next time step call to integrate. This represents 
the minimum amount of communication necessary between tasks in order to ensure the code 
runs like the serial version. In order to achieve a good balance of computational workload 
between tasks, the code dynamically adjusts how many radial bins are devoted to each task 
based on how much time the previous time step took over the various radial ranges. This 
approach to optimization is based on the assumption that the system will not change too 
drastically each At. Given that the time steps we use are relatively small, the radial ranges 
change only gradually, and thus the resulting computations remain quite similar. 

The amount of speed up in code execution achieved through our simplest parallelization 
approach compared to the original BT code of Durisen and colleagues scales as the number 
of processors (~ n p ) on a single node since the number of operations involved in computing 
the loss and gain integrals scales as the number of radial bins ( N ). Typically a single node 
(e.g., a quad core CPU) will have n p = 8 processors, which is what we have used for the bulk 
of the simulations for this paper. We have done some cases using a CPU with n p = 12. Thus, 
using the minimal communication parallelization, we achieve a speedup of a factor of ~ 10 
compared to the serial BT code. We have developed a fully parallelized version of the BT 
code (maximum amount of communication), which will likely be useful when expanding the 
size of the radial grid (i.e., global simulations with many radial bins) and requiring multiple 
nodes to do so. Multiple nodes naturally introduces another layer of communication to 
contend with, communication between nodes in addition to tasks within a node. However, 
for the demonstrative simulations we present in this paper, the minimum communication 
BT code is more than sufficient. 
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Figure 1: The optical depth profile from the UVIS Cassini ct-Arae occultation (Colwell 
et al, 2009) as a function of radius from Saturn in planetary radii. Top panel: Much of 
the structure in the C ring has long remained a mystery; however, there is evidence to 
suggest ballistic transport is at work, particularly with the inner B ring edge sharpness 
and its associated ramp, and the outer (and inner) edges of plateaus. The radial drift due 
to ballistic transport can lead to a pile up of material at plateau edges, which may be 
characterized by the sloped tops of the plateaus. The reversal of the slopes of plateau peaks 
about the Maxwell ringlet at ~ 1.45 Rg may be indicative of impacts that change in nature 
to fragmenting from cratering across the region. Bottom panel: The Cassini division and 
inner A ring displays some of the same properties as the C/inner B ring. A similar ramp 
connects to a sharp inner edge. Even though the radial scales in the plots are different, the 
ramp width is roughly the same in both cases. 

Figure 2: Reproduction of Fig. 10a from Durisen et al. (1992) where Y — 3 x 10 5 and 
Xb = 2 x 10~ 4 . The run spans 99 to with t(R) plotted at selected times (black curves). Also 
plotted are the associated fractional masses of pollutant (red curves) for these structural 
runs using / ext = 0.5 and rj = 0.12. Because the opacity is assumed to be constant, there is 
not much contrast between the C ring and inner B ring in composition at longer times. Our 
structural results are identical to Durisen et al. (1992); the edge is sharpening and a ramp 
is emerging (Sec. 3.1). 

Figure 3: Reproduction of Fig. 10b from Durisen et al. (1992) which details the ramp and 
B-C boundary region further inwards into the C ring at different times (black curves). The 
compositional evolution is shown in red. Although there is more contrast in composition for 
the longest time here, the variation is due to brighter material spilling over from the inner B 
ring edge into the C ring. The net contrast across the 1.49 — 1.54 R^ region remains roughly 
constant. The undulations seen may come as a result of the BTI. 

Figure 4: Reproduction of Fig. 15 from Durisen et al. (1992) which is a compilation of three 
different simulations for Y — 3 x 10 5 (dashed curve), 10 6 (dotted curve) and 2 x 10 6 (solid 
curve) over a simulation time of 104 to- The models are plotted against the Voyager PPS 
imaging optical depth. The larger yield models use a smaller Xb = 10” 4 and have higher 
resolution. Our structural results are essentially identical to Durisen et al. (1992). The 
compositional models are shown in red. In this case, the strongest compositional gradient is 
seen in the highest yield case because of the much shorter absolute timescale involved (also 
discussed by CE98). 

Figure 5: A complementary suite of simulations to Fig. 4 in which we have employed the 
variable opacity model (Eq. [22]). This model gives a slightly larger mass in the inner B ring, 
and more mass in the plateaus, but five times less mass, overall, in the C ring compared 
to the constant opacity model. This is quite evident in the compositional evolution (red 
curves) where now much larger absolute contrasts in fractional mass are seen compared to 
the constant opacity model (the scale is the same as Fig. 4). The major structural difference 
after 104 t^ is in the inner B ring where, in particular, the undulatory structure for the 
2 x 10 6 case is quite suppressed. 
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Figure 6: A plot of the time rate of change of optical depth as a function of ring radius 
for the Y — 2 x 10 6 model in Fig. 4. These plots break down the individual contributions 
to structural change from purely BT effects such as mass transfer and advection (of BT 
only) in the top panel, and the effect of viscosity on the total in the bottom panel. The 
first undulation (the steady-state edge) is visible around R ~ 1.525 R 5 . Outside of this 
location, the growth of quasi-steady-state undulations in the inner B ring is due to a lack of 
cancellation between BT and viscosity (the net effect is plotted in red). This actually leads 
to uninhibited growth in amplitude over much longer timescales. On the other hand, the 
strong variations inside the inner B ring edge, which are due almost entirely to BT, tend 
to vary in location and thus the effects tend to average out over time so that similar large 
amplitude growth of structure is not seen. However, these variations do lead to structure at 
later times (e.g., see Fig. 8 ). 

Figure 7: A plot of the time rate of change of optical depth as a function of ring radius 
for the Y = 2 x 10 6 model in Fig. 5 which uses the variable opacity model. In this case, 
the lack of undulations in the inner B ring is explained by both the smaller amplitude of the 
effects of BT and viscosity, and a higher level of cancellation between them (red curve). A 
combination of slightly more mass in the inner B ring, as well as much less mass in the C 
ring, and the variable k combine to prevent large amplitude growth. 

Figure 8: A long term simulation of the structural evolution of the inner B ring edge for 
Y = 1.5 x 10 6 , Xb = 10 ” 4 and f v = 1 for a constant opacity model. The total simulation 
time is 800 t <3 which corresponds to an absolute timescale of ~ 7 x 10 6 years for this set 
of parameters. The simulation clearly shows how BT can maintain a sharp edge for long 
periods. A ramp forms, but develops a great deal of structure which is not present in the 
observations (dot-dashed curve). Also shown is the 800 tc profile (red curve) for the variable 
opacity model. See Sec. 3.2 for more details. 

Figure 9: Suite of models for Y — 3 x 10 5 , Xb = 10 ~ 4 and f u = 0.125 that employ different 
steepness in the slope of the ejecta velocity distribution (Eq. [25], n — 3 is our adopted 
fiducial value) plotted against the a-Arae occupation r (dotted curve). The variation leads 
to two effects. First, lower n sharpens the inner edge much more effectively than our fiducial 
value (for the same parameter choices) while seemingly producing a large, featureless ramp 
(although with less steepness). Second, larger n produce “notch” like features in the edge 
not unlike that seen in the real optical depth profile. 

Figure 10: Two hypothetical massive B ring models in which both structure (black curves) 
and pollution transport (red curves) evolve over time (100 tc). The variable opacity model is 
used along with a modified kinematic viscosity in which the magnitude of the viscosity at the 
inner and outer B ring edges is chosen to maintain the edges close to their initial sharpness. 
The dashed curves correspond to a model in which the B ring core (highest optical regions 
between ~ 1.65 — 1.85 R 5 ) has a ~ 400 g cm -2 ; the dotted curves have a core a = 1000 
g cm -2 . A significant difference in the core composition is seen between the two models. 
The traditional “Mimas” mass model is also shown (solid curves) for comparison. The initial 
condition for the Mimas mass model only is given by the blue curve. For more details, see 
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Sec. 3.2. 


Figure 11 : Compilation of combined structural and compositional models with Y = 10 4 
(solid curves), 10 5 (dotted curves) and 10 6 (dashed curves) using the u-Arae occupation r 
as the initial condition (blue curve). The simulations are run over 2 t^ in the active radial 
range of 1.4 — 1.53 Rg and assume Xb = 10 -4 , a Wisdom and Tremaine kinematic viscosity 
with f v = 1 and a variable opacity model. See Sec. 3.3 for explanation. 

Figure 12: Evolution of plateau features centered at 1.465 R# using a constant Y / f u = 10 5 
over 2 tc that demonstrates the role of the kinematic viscosity in maitaining ring structure. 
For the yields of 10 4 (black curve) and 10 6 (blue curve), the dominant particle size f is defined 
to be 20.4 and 204 cm, respectively. Despite the large variation in timescale associated with 
the different yields, the profiles look similar due to the degree of scaling that exists between f„ 
and Y (see text). Also due to the different timescales, there has been some radial registration 
of the profiles to line them up for easier comparison. 

Figure 13: Both the structural and compositional evolution of the same plateaus as Fig. 
12 for Y = 10 5 and f u = 0.3 ( Xb = 10~ 4 , black curve). Here we vary the lower bound Xb of 
an n — 3 ejecta velocity distribution. The blue curve corresponds to Xb = 5 x 10 -5 (u e j ~ 1 
m s~ 4 ), the red curve to Xb = 2 x ICE 5 (~ 0.4 m s” 1 ) and the green curve to Xb = 2 x 10 -4 
(~ 4 m s _1 ). The effect of lowering Xb is somewhat similar to increasing the viscosity (BT is 
weaker because of a shorter throw distance). In this case, Xb = 5 x ICE 5 appears to provide 
a better fit to the plateaus’ outer edges. The largest Xb demonstrates how the larger throw 
distances (and thus angular momentum disparities) lead to more pileup of material along 
with more radial drift. In the lower panel, the compositional profiles with smaller Xb show the 
most variation because the bulk of material does not travel as far from its point of ejection. 

Figure 14 : Evolution over 2 R; of plateau features centered at 1.495 Rs with Y = 10 6 , and 
Xb = ICE 4 for different viscosity models (Sec. 2.2.7). The black curve corresponds to a purely 
kinematic viscosity model, whereas the red and blue curves correspond to combined models 
with different ring particle densities that include a viscosity due to wake formation. In all 
cases, f u — 1 for the kinematic viscosity. The condition for wake formation is marginally 
satisfied in the C ring plateaus. Compared to the purely kinematic viscosity model, the 
higher viscosity in the plateaus is very effective at suppressing the large amplitude growth 
of structure there despite such a high yield. 

Figure A-l: Illustration that demonstrates the relationship between the “radiative transfer” 
angles ( 6 , f) and the cone and clock angles of emergent ejecta, ( a , (3 ). The angles are defined 
in a frame rotating with the ring material where a is measured from the orbital velocity 
vector Vk, and /3 from the ring normal. Reproduced from CD90. 

Figure A-2: An example of actual model calculations of the ejecta distribution function 
& at R = 1.8 Rs, which is plotted for two different optical depths as a function of the 
cone and clock angles a and f3 (in radians) measured from the ring particle velocity vector. 
These models are averaged over ring orbit longitude. The prograde bias of the distribution 
is illustrated by the maximum occurring for a < 1.5. The parameterized fit we use in the 
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BT code (Eq. [A-l] ) matches the actual calculations within 10%. Reproduced from CD90. 

Figure C-l: Plot of the magnitude of the theoretical viscosity as a function of radius in the 
rings for a purely kinematic viscosity model (bold curve, Eq. [24]) which is r-dependent, and 
the viscosity model due to wakes (when greater, lighter curve, Eq. [C-l]) which depends on 
cr. The UVIS Cassini a-Arae occupation data defines r, and a variable opacity is used. The 
density of ring particles for the wake viscosity is assumed to be 0.5 g cm -3 , while f = 64.6 
cm for the kinematic viscosity. The condition for wake formation can be marginally satisfied 
in the C ring plateaus, especially if the mean density of particles is higher. Outside of the 
plateaus, the viscosity is purely kinematic. Some evidence exists for wakes in the inner B 
ring, which would mean much larger viscosities than given by a purely kinematic model. 
This could pose a problem for maintaining a sharp edge by BT alone, but the situation is 
complicated (see Sec. 2.2.7). 
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Table 1: Model Parameters 


Parameter 

Definition 

sY 

Probability of meteoroid impact for given r 

B 

Elevation angle 

D 

Edge width parameter for inner B ring edge 
models 

f 

Ejecta speed distribution 

fe 

Total mass fraction of extrinsic pollutant 

/ext 

Mass fraction of extrinsic, absorbing material 

Fg 

Gravitational focusing factor 

u 

Kinematic viscosity scaling factor 

h, h c 

Specific angular momentum per unit mass 

V A 

1 m,hj 1 l m,h 

Direct gains and losses in mass and angular 
momentum 

r r • 

L m,e; L m,i 

Fractional mass gain integrals for extrinsic and 
intrinsic materials 

V 

Retention efficiency of absorbing properties fol- 
lowing impact 

K 

Opacity = r/cr 

V 

Ring viscosity 

N 

Number of Lagrangian (annuli) bins 

n 

Power law exponent for / 

& 

Ejecta absorption probability 

r 

Ring particle radius 

n,r 2 

Upper and lower bound of particle size 
distribution 

f 

Mean particle radius above which 1/3 of the op- 
tical depth is accounted for 

R 

Radial distance from the central planet 

P 

Ring particle density 

a 

Surface mass density 

d e j 

Ejecta mass flux 


Two sided mass flux impacting within the rings 

doo 

Unfocussed, one-sided, interplanetary mass flux 

r 

Normal optical depth 

tc 

Gross erosion time 

vr 

Total radial drift velocity 

X 

Ratio of ejecta velocity to orbital velocity 

Xb 1 Xt 

Lower and upper bounds of ejecta distribution 

Y 

Impact ejecta yield 


Distribution of ejecta with angle and optical depth 
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Figure 1: Estrada et al. 
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Figure 2: Estrada et al. 
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Figure 3: Estrada et al. 
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Figure 4: Estrada et al. 
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Figure 5: Estrada et al. 
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Figure 6: Estrada et al. 
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Figure 7: Estrada et al. 
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Figure 8: Estrada et al. 
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Figure 9: Estrada et al. 
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Figure 10: Estrada et al. 
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Figure 11: Estrada et al. 


62 


fractional mass 


optical depth 



ring radius (R s ) 


Figure 12: Estrada et al. 
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Figure 13: Estrada et al. 
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Figure 14: Estrada et al. 
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